% Harj. 6 teht. 5  (ratk6teht5s.m)
function ratk6teht5s

opt=odeset('Jacobian',@jakobiaani);
opt=odeset(opt,'abstol',1e-5,'reltol',1e-5);
opt=odeset(opt,'Stats','on');
['abstol=',num2str(opt.AbsTol)]
['reltol=',num2str(opt.RelTol)]


tspan=[0 3]; y0=[1;0;0];
%global kutsut
%kutsut=0
%[t45,y45]=ode45(@kemikaalit,tspan,y0,opt);   % Kestää jonkin aikaa
%subplot(121), plot(t45,y45(:,2),'r')
%aks=axis; aks(1)-0.2; axis(aks);
%xlabel('t'), ylabel('y_2(t)'), title('ode45','FontSize',14)


[ts,ys]=ode15s(@kemikaalit,tspan,y0,opt);


subplot(3,1,1)
plot(ts,ys(:,1),'b')
xlabel('t'), ylabel('y_1(t)'), title('ode115s','FontSize',14)

subplot(3,1,2)
plot(ts,ys(:,2),'b')
xlabel('t'), ylabel('y_2(t)'), title('ode115s','FontSize',14)

subplot(3,1,3)
plot(ts,ys(:,3),'b')
xlabel('t'), ylabel('y_3(t)'), title('ode115s','FontSize',14)

figure(2)
plot3(ys(:,1),ys(:,2),ys(:,3))
grid on

function yp=kemikaalit(t,y)
y1=y(1);y2=y(2);y3=y(3);
yp=[-0.04*y1+10^4*y2*y3;
    0.04*y1-10^4*y2*y3-3*10^7*y2;
    3*10^7*y2^2];



function z=jakobiaani(t,y)
y1=y(1);y2=y(2);y3=y(3);
z=[-.4e-1, 10000*y3, 10000*y2; 
    .4e-1, -10000*y3-30000000, -10000*y2; 
    0, 60000000*y2, 0];

