clear;
xset("window",0);
function z=f(y); z=-y; endfunction;
function z=ff(t,y); z=f(y); endfunction
function z=fff(t,v); z=[1,f(v(2))]'; endfunction;
fchamp(fff,0,0:0.2:2,0:0.4:4);

Tmax=2;y0=4;
N=100;smallstep=Tmax/N;
M=10;largestep=Tmax/M;
t=0:smallstep:Tmax;
y=ode(y0,0,t,ff);
plot(t,y);

//Euler(t0,y0,step) is an approxroximation
//of y(t0+step) pour y(t0)=y0.
function y=Euler(t0,y0,step);
  y=y0+step*ff(t0,y0);
endfunction;

plot([0,largestep],[y0, Euler(0,y0,largestep)],'r-o');
y1=Euler(0,y0,largestep)
plot([largestep,2*largestep],[y1,Euler(largestep,y1,largestep)],'r-o');

step=largestep;
tt=0:step:Tmax;
yy=zeros(1+Tmax/step);yy(1)=y0;
for n=1:1:Tmax/step;
  yy(n+1)=Euler(tt(n),yy(n),step);
end;
plot(tt,yy,'g->'); // "g" for "green"
disp(y(N+1)-yy(1+Tmax/step),'difference=',step,'lstep=');

// now using a small step
step=smallstep;
tt=0:step:Tmax;
yy=zeros(1+Tmax/step);yy(1)=y0;
for n=1:1:Tmax/step;
  yy(n+1)=Euler(tt(n),yy(n),step);
end;
plot(tt,yy,'r-.');
disp(y(N+1)-yy(1+Tmax/step),'difference=',step,'step=');

///////////////////////// now with -y^2
xset("window",1);
function y=f(y); y=-y.*y; endfunction;
function y=ff(t,y); y=f(y); endfunction;

Tmax=2;y0=4;
N=100;smallstep=Tmax/N;
M=10;largestep=Tmax/M;

t=0:smallstep:Tmax;
y=ode(y0,0,t,ff);
plot(t,y);

///@
function y=Euler(t0,y0,step);
  y=y0+step*ff(t0,y0);
endfunction;

plot([0,largestep],[y0, Euler(0,y0,largestep)],'r-o');

step=largestep;
tt=0:step:Tmax;
yy=zeros(1+Tmax/step);yy(1)=y0;
for n=1:1:Tmax/step;
  yy(n+1)=Euler(tt(n),yy(n),step);
end;
plot(tt,yy,'g->');
disp(y(N+1)-yy(1+Tmax/step),'difference=',step,'step=');

// now with a small step
step=smallstep;
tt=0:step:Tmax;
yy=zeros(1+Tmax/step);yy(1)=y0;
for n=1:1:Tmax/step;
  yy(n+1)=Euler(tt(n),yy(n),step);
end;
plot(tt,yy,'r-');
disp(y(N+1)-yy(1+Tmax/step),'difference=',step,'step=');

///////////////////////// on recommence avec -y^1/2
xset("window",2);
function y=f(y); y=-sqrt(abs(y)); endfunction;
function y=ff(t,y); y=f(y); endfunction

Tmax=2;y0=4;
N=100;psmallstep=Tmax/N;
M=10;largestep=Tmax/M;

t=0:smallstep:Tmax;
y=ode(y0,0,t,ff);
plot(t,y);


function y=Euler(t0,y0,step);
  y=y0+step*ff(t0,y0);
endfunction;

plot([0,largestep],[y0, Euler(0,y0,largestep)],'r-o');

step=largestep;
tt=0:step:Tmax;
yy=zeros(1+Tmax/step);yy(1)=y0;
for n=1:1:Tmax/step;
  yy(n+1)=Euler(tt(n),yy(n),step);
end;
plot(tt,yy,'g->');
disp(y(N+1)-yy(1+Tmax/step),'difference=',step,'step=');

// now with a smallstep
step=smallstep;
tt=0:step:Tmax;
yy=zeros(1+Tmax/step);yy(1)=y0;
for n=1:1:Tmax/step;
  yy(n+1)=Euler(tt(n),yy(n),step);
end;
plot(tt,yy,'r-');
disp(y(N+1)-yy(1+Tmax/step),'difference=',step,'step=');

/////////// again y'=-y with RK2 method
xset("window",3);
// fchamp()
xr=0:0.2:2;yr=0:0.4:4;
//function g=ggg(t,v); 
//  g(1)=1;g(2)=-v(2); endfunction;
//  fchamp(ggg,0,xr,yr);
function ff=f(y); ff=-y; endfunction;
function ff=ff(t,y); ff=f(y); endfunction
function ff=fff(t,v); ff=[1,f(v(2))]';// v is 2-dimensional
endfunction;
xr=0:0.2:2;yr=0:0.4:4;
fchamp(fff,0,xr,yr);

Tmax=2;y0=4;
N=100;smallstep=Tmax/N;
M=10;largestep=Tmax/M;

t=0:smallstep:Tmax;
y=ode(y0,0,t,ff);
plot(t,y);

//RK(t0,y0,step) is an approximation of y(t0+pas) pfor y(t0)=y0.
function y=RK(t0,y0,step);
//  local k1,k2;
  k1=step*ff(t0,y0);
  k2=step*ff(t0+step/2,y0+k1/2);
  y=y0+k2;
endfunction;

plot([0,largestep],[y0, RK(0,y0,largestep)],'r-o');

step=largestep;
tt=0:step:Tmax;
yy=zeros(1+Tmax/step);yy(1)=y0;
for n=1:1:Tmax/step;
  yy(n+1)=RK(tt(n),yy(n),step);
end;
plot(tt,yy,'g->');
disp(y(N+1)-yy(1+Tmax/step),'difference=',step,'step=');

// and now with a small step
step=smallstep;
tt=0:step:Tmax;
yy=zeros(1+Tmax/step);yy(1)=y0;
for n=1:1:Tmax/step;
  yy(n+1)=RK(tt(n),yy(n),step);
end;
plot(tt,yy,'r-.');
disp(y(N+1)-yy(1+Tmax/step),'difference=',step,'step=');
disp(4*exp(-2)-y(N+1),'difference=4exp(-2)-ode(2)');  
disp(y(N+1)-yy(1+Tmax/step),'difference=',step,'step=');//step=Tmax/M;
