clear

sysName = 'Pendulum';
      S = eval(['sys',sysName,'([])']);
    dim = S.dimQ;
      T = 1000;
options = optimset(optimset('fsolve'),'MaxIter',1000,'TolFun',1e-15, ...
                   'TolX',1e-15);

  mtd = 'euler';
order = 1;
 hmin = eval(['1e-',int2str(order)]);
 Nmax = int32(T/hmin);
 file = ['data',sysName,'/',mtd,'_logN',int2str(order),'_T',int2str(T), ...
         '.mat'];
if exist(file,'file')
    load(file)
else
     f = @(t,z,S) [z(dim+1:2*dim); S.ELeq(z(1:dim),z(dim+1:2*dim),S)];
    zs = ode_rk2(f,0:hmin:T,[S.q0;S.v0],S);
    save(file,'T','order','hmin','Nmax','zs')
end

for s=1:5

order = 1;
    h = eval(['1e-',int2str(order)]);
    N = int32(T/h);
%    s = 4;
 file = ['data',sysName,'/vi_s',int2str(s),'_logN',int2str(order), ...
         '_T',int2str(T),'.mat'];
if exist(file,'file')
    load(file)
else
    [q,p] = varInt(S,S.q0,S.p0,T,N,s,options);
    save(file,'T','order','h','N','s','q','p')
end
q = reshape(q(:,1,:),[N+1,1]);

format shorte;
disp(max(abs(zs(1,1:Nmax/N:end)-q')))

end