function [dy]=oscillateur(t,y),
	dy(1)=y(2),
	dy(2)=-4*%pi^2*y(1),
endfunction 
			
function [y]=pasEuler(systeme,ne,h,t,yn),
			dy = systeme(t,yn);
			for e=1:ne,
				y(e) = yn(e)+h*dy(e)
			end;
endfunction
			
			
function [mt,my]=euler(yi,T,N,systeme,ne),
			h = T/N;
			my = zeros(ne,N+1);
			mt = zeros(1,N+1);
			my(:,1)=yi;
			mt(1)=0;
			for k=1:N,
				my(:,k+1)=pasEuler(systeme,ne,h,mt(k),my(:,k));
				mt(k+1) = mt(k)+h;
			end;
endfunction
			
yi=[1;0];
ne=2;
N=1000;
T=5;
[mt,my] = euler(yi,T,N,oscillateur,ne);
			
figA=scf();
plot2d(mt,my(1,:));
xtitle('','t','y(1)');
			
figB=scf();
plot2d(mt,abs(my(1,:)-cos(2*%pi*mt)));
xtitle('','t','erreur');
			

N=N*2;
[mt2,my2] = euler(yi,T,N,oscillateur,ne);
N=N*2;
[mt3,my3] = euler(yi,T,N,oscillateur,ne);
			
figC=scf();
plot2d(mt,abs(my(1,:)-cos(2*%pi*mt)),style=4);
plot2d(mt2,abs(my2(1,:)-cos(2*%pi*mt2)),style=5);
plot2d(mt3,abs(my3(1,:)-cos(2*%pi*mt3)),style=6);
xtitle('','t','erreur');
			
function my=eulerEnd(yi,T,N,systeme,ne),
			y = yi;
			h = T/N;
			t = 0;
			for k=1:N,
				my = pasEuler(systeme,ne,h,t,y);
				y = my;
				t = t+h;
			end;
endfunction			
			
N=logspace(2,5,30);
e=N;
h=N;
k=size(N); k=k(2);
T=4.5;
for i=1:k,
	my = eulerEnd(yi,T,N(i),oscillateur,ne);
	e(i) = abs(my(1)-cos(2*%pi*T));
	h(i) = T/N(i);
end;
			
figD=scf();
plot2d(h,e,style=-5,logflag='ll');
xtitle('','h','erreur');
			
T=10;
N=10000;
[mt,my] = euler(yi,T,N,oscillateur,ne);
			
figE=scf();
plot2d(my(1,:),my(2,:)/(2*%pi));
xtitle('','y(1)','y(2)');
			
M=[0,1;-4*%pi^2,0];
vp = spec(M);
			
M=[0,1;-4*%pi^2,-0.1];
vp = spec(M);
			
n=20;
h=logspace(-5,-2,n);
m = h;
for i=1:n,
	m(i)=(1-0.05*h(i))^2+4*%pi^2*h(i)^2;
end;
			
figF=scf();
plot2d(h,m,logflag='ln');
xtitle('','h','module')
			
function [dy]=amorti(t,y),
	dy(1)=y(2),
	dy(2)=-4*%pi^2*y(1)-0.1*y(2),
endfunction
			

T=30;
h=2e-4;
N=int(T/h);
ne=2;
yi=[1;0];
[mt1,my1] = euler(yi,T,N,amorti,ne);
			
figG=scf();
plot2d(my1(1,:),my1(2,:)/(2*%pi));
xtitle('','y(1)','y(2)');
			

h=4e-3;
N=int(T/h);
ne=2;
yi=[1;0];
[mt2,my2] = euler(yi,T,N,amorti,ne);
			
figH=scf();
plot2d(my2(1,:),my2(2,:)/(2*%pi));
xtitle('','y(1)','y(2)');
			
s=size(mt1); s=s(2);
e=zeros(s);
for i=1:s,
	e(i) = abs(my1(1,i)-exp(-0.05*mt1(i))*(cos(2*%pi*mt1(i))+0.05/(2*%pi)*sin(2*%pi*mt1(i))));
end;
		
figJ=scf();
plot2d(mt1,e);
xtitle('','t','erreur');
		
