Initialize
x_old = 43100, y_old = 0 and z_old = 0.0033
Compute x_new by solving the nonlinear system of equations
(x_new-x_old)/dt = lambda - rho*x_new - beta*x_new*z_new
(y_new-y_old)/dt = beta*x_new*z_new - delta*y_new
(z_new-z_old)/dt = p*y_new - c*z_new
by fixed-point iteration or with MATLAB's fsolve, e.g.
This gives you the solution for your system at time t=dt.
Set
x_old = x_new, y_old = y_new and z_old = z_new
and solve the above system again for x_new, y_new and z_new.
This gives you the solution at time t=2*dt.
Continue until you reach t=tfinal.
Best wishes
Torsten.
