Hello, I have to write a code in C which depicts Philip Munz's "When Zombies Attack" population model. I have attached the link to this article if you want to read it but the basic idea is that I have to solve a system of 3 first order differential equations and a system of 4 first order differential equations in C which are:

3 first order system of ODE's (SZR model):

S' = SZ

Z' = SZ + R SZ

R' = S + SZ R .

4 first order system of ODE's (SIZR model)

S' = SZ

I' = SZ I I

Z' = I + R SZ

R' = S + I + SZ R

Where ...

S is Susceptible Humans

Z is Zombies

R is Removed individuals (people who have died by zombies or natural causes)

I is the infected class of people who remain this way for a period of time before they turn into a zombie.

I have started a code in C and need help to finish it. My task is to run Munz's basic (SZR) and latent infection (SIZR) models using the rk4 ODE solving method code which I have partially written in Matlab and got stuck. I have provided all the codes that I have started to help explain what I need to accomplish. Please help me understand how to finish my code in C. I do not need it finished in Matlab I only used this to help me understand.

Matlab Code

%% "MAIN" h = 0.01; t = [0:h:20]; y = zeros(length(t),3);

y(:,1) = 500; %S y(:,2) = 0; %Z y(:,3) = 0; %R

onesix = 1/6;

% theta -> y(1); % omega -> y(2);

for i = 2: 1: length(t) y(i,:) = myrk4(t(i), y(i-1,:), h); end

plot(t,y(:,1)); hold on; plot(t,y(:,2));

if exist('Soln.txt', 'file') fd = fopen('Soln.txt', 'r'); ODEsoln = textscan (fd, '%f %f %f', length(t)); fclose(fd); tC = ODEsoln{1}; xC = ODEsoln{2}; plot(tC, xC); else disp('No C results for comparison'); end

%% RK4 Step: myrk4 function [ yout ] = myrk4(t, y, h) hh = h*0.5; onesix = 0.1666666666666667;

dy1 = h * CMO_der(t, y); dy2 = h * CMO_der(t+hh, y+0.5*dy1); dy3 = h * CMO_der(t+hh, y+0.5*dy2); dy4 = h * CMO_der(t+h, y+dy3);

yout = y + onesix*(dy1+2*dy2+2*dy3+dy4); end

%% "Derivs" Function, specific to pendulum function [ dy ] = CMO_der(t, y) a = 0.005; b = 0.0095; ze = 0.0001; d = 0.0001; dy(1) = b*y(1)*y(2); dy(2) = b*y(1)*y(2) + ze*y(3) - a*y(1)*y(2); dy(3) = d*y(1) + a*y(1)*y(2) - ze*y(3); end

C Code Stub:

#include #include #include

#define STEP_SIZE 0.01 #define END_TIME 10.0 #define NUM_EQNS 2

void rk4_step (int, double, double, double*, double*); void derivs (int, double, double*, double*);

int main() { double **y, *t; int i, j, nsteps; FILE* TxtOut;

nsteps = END_TIME/STEP_SIZE+1;

/* Allocate t[nsteps], y[nsteps][NUM_EQNS] */

/* Initial Condition: Susceptibles = 500 */ y[0][0] = 500.; t[0] = 0.;

/* Start Main Integration Loop */ for (i=1; i

/* File I/O */ if ( (TxtOut = fopen("Soln.txt", "w")) == NULL ) { fprintf(stderr, "Cannot open matrix text file OUTPUT... exiting "); exit(-1); } for (i = 0; i < nsteps; i++) { fprintf (TxtOut, "%18.14e ", t[i]); for (j = 0; j < NUM_EQNS; j++) { fprintf (TxtOut, "%18.14e ", y[i][j]); } fprintf (TxtOut, " "); }

if (fclose(TxtOut) == EOF) { fprintf(stderr, "Cannot close matrix text file INPUT... exiting "); exit(-1); }

return(0); }

void rk4_step (){ }

void derivs ( ) {

}

rk4 code in C:

#include #include #include

double rk4(double(*f)(double, double), double h, double t, double y) { double k1 = h * f(t, y), k2 = h * f(t + h / 2, y + k1 / 2), k3 = h * f(t + h / 2, y + k2 / 2), k4 = h * f(t + h, y + k3);

return y + (k1 + 2 * k2 + 2 * k3 + k4) / 6; }

double rate(double t, double y) { return t * sqrt(y); }

int main(void) { double *yout, t, y2; double t0 = 0, t1 = 10, h = .1; int i, n = 1 + (t1 - t0)/h; yout = malloc(sizeof(double) * n);

for (yout[0] = 1, i = 1; i < n; i++) yout[i] = rk4(rate, h, t0 + h * (i - 1), yout[i-1]);

printf("x\tyout\trel. err. ------------ "); for (i = 0; i < n; i += 10) { t = t0 + h * i; y2 = pow(t * t / 4 + 1, 2); printf("%g\t%g\t%g ", t, yout[i], yout[i]/y2 - 1); }

return 0; }