1. The problem statement, all variables and given/known data Hi! I'm to build a model which solves the basic dynamic eqns in the so called SIR model in epidemology. I'm using MatLab and RK4. There are tree populations in this model: Succeptible S, infected I and recovered R. 2. Relevant equations The eqations are dS/dt = -beta * S*I/N dI/dt = beta * S*I/N - alpha*I dR/dt = alpha*I where alpha and beta are rates of recovery and infection (constants). Vaccination will be introduced aswell, but I have set that to 0 for now. 3. The attempt at a solution My base .m-file is this: Code (Text): clear all close all %Step size and number of iterations are defined dt=1; N=1000; %days P = 1000000; %The initial conditions are entered as the first row in the array z. z(1)=P-1; % = S z(2)=1; % = I z(3)=1; % = R t(1)=0; %RK4 on function for n=1:(N) k1=P1f(t(n),z(n,:)); k2=P1f(t(n)+dt/2,z(n,:)+1/2*k1*dt); k3=P1f(t(n)+dt/2,z(n,:)+1/2*k2*dt); k4=P1f(t(n)+dt,z(n,:)+k3*dt); z(n+1,:)=z(n,:)+1/6*(k1+2*k2+2*k3+k4)*dt; t(n+1)=t(n)+dt; end %The resluts are plotted subplot(1,2,2) plot(t,z(:,2),'r') axis([0 1000 0 100000]) xlabel(''),ylabel('') subplot(1,2,1) plot(t,z(:,3)) axis([0 1000 0 100000]) xlabel(''),ylabel('') with the function P1f.m: Code (Text): function f=P1f(t,z) %Def constants alpha = 0.25; sigma = 50; T = 0.1; beta = sigma*T; gamma= 0; f(1) = -beta*z(1)*z(2)/z(3) - gamma*z(1); f(2) = beta*z(1)*z(2)/z(3) - alpha*z(2); f(3) = gamma*z(1) + alpha*z(2); end The constants k in the RK4 routine all become NaN, and I don't understand why. This happenes during the first iteration. Most of my problem is that I don't fully understand how the 'function' definition in matlab works. As you can see, I only define the RHS of the original eqations for the first 3 elements in z. How does the function calculate the values past that first row? I know the method is sound because I've used it before, last year, and I've forgotten how it works. Very grateful for any help, Eric EDIT: Major mistake was dividing by R instead of N = R + I + S in the right hand sides. Seems to be working now, as far as I can see.