f = @(t,y,s) (0.0012*((1-y)/(1.66*10^(6)))+s-0.0012*0.02*y); %y=Nddotexi,vmed=0.02,fdot=3*10^(-18),r=0.0012,k=Nsurv
a = 0; %input('initial ponit, a: ');
b = 100; %input('end point, b: ');
n = 10; %input('intervals, n: ');
alpha = 0; %input('initial condition, alpha: ');
h = (b-a)/n;
t=[a zeros(1,n+1)];
y=[alpha zeros(1,n+1)];
s=[0 zeros(1,n+1)];
for i = 1:n+1
t(i+1)=t(i)+h;
yprime=y(i)+h*f(t(i),y(i),s(i)); %slope
y(i+1)=y(i)+h*f(t(i+1),yprime,s(i+1));
fprintf('%.4f %.4f\n', t(i), y(i));
end
figure(1)
plot(t,y,'r-o');
xlabel('t values'); ylabel('y values');
grid on;
legend('y')
