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
			
function [J]=jacobienOscillateur(t,y),
			J(1,1) = 0; J(1,2)=1;
			J(2,1)=-4*%pi^2; J(2,2)=0;
endfunction			
			
function [y]=pasEulerImplicite(systeme,jacobien,Id,P,h,t,yn,tol),
			stop = 0; it=0; maxit=100;
			y = yn;
			t1 = t+h;
			while stop==0,
				it = it+1;
				dy = systeme(t1,y);
				G = yn+h*dy-y;
				[deltaY,ker] = linsolve(h*jacobien(t1,y)-Id,G);
				y = y+deltaY;
				if it>maxit then stop=1; end;
				stop = 1;
				for k=1:P,
					if abs(deltaY(k))>tol then stop=0; end;
				end;
			end;
endfunction
			
function [mt,my]=eulerImplicite(yi,T,N,systeme,jacobien,P,tol),
			Id = eye(P,P);
			h = T/N;
			my = zeros(P,N+1); 
			mt = zeros(1,N+1);
			my(:,1)=yi;
			mt(1)=0;
			for k=1:N,
				my(:,k+1)=pasEulerImplicite(systeme,jacobien,Id,P,h,mt(k),my(:,k),tol);
				mt(k+1) = mt(k)+h;
			end;
endfunction
			
yi=[1;0];
P=2;
N=5000;
T=20;
tol=1e-5;
[mt,my] = eulerImplicite(yi,T,N,oscillateur,jacobienOscillateur,P,tol);
			
figA=scf();
plot2d(mt,my(1,:));
xtitle('','t','y(1)');
			
yi=[1;0];
N=10000;
[mt,my] = eulerImplicite(yi,T,N,oscillateur,jacobienOscillateur,P,tol);
			
figB=scf();
plot2d(mt,my(1,:));
xtitle('','t','y(1)');
			
function [dy]=amorti(t,y),
	dy(1)=y(2),
	dy(2)=-4*%pi^2*y(1)-0.1*y(2),
endfunction

function [J]=jacobienAmorti(t,y),
	J(1,1)=0; J(1,2)=1;
	J(2,1)=-4*%pi^2; J(2,2)=-0.1;
endfunction

P=2;
N=10000;
T=40;
tol=1e-5;
[mt,my] = eulerImplicite(yi,T,N,amorti,jacobienAmorti,P,tol);
function y=f(t),
	y=exp(-0.05*t)*(cos(2*%pi*t)+0.05/(2*%pi)*sin(2*%pi*t))
endfunction
[mt1,my1] = euler(yi,T,N,amorti,P);
			
figC=scf();
plot2d(mt,my(1,:),style=color('red'));
plot2d(mt1,my1(1,:),style=color('blue'))
fplot2d(mt,f)
xtitle('','t','y(1)');
			
