Runge kutta asignment Matlab

Question

I need help with a home asignment I have for a matlab course. The aim is to solve P and Z for the two functions dP/dt and dZ/dt. The asignment is to a) solve for P and Z. b) Plot the error function for P and Z by calculating |P(n)h/2 - P(n)h| where P(n)h is P with regular steplength, P(n)h/2 is steplength halved. dP/dt = rP (1 - P2/α2 + P2) - RmZ dZ/dt = γRmZ P/(K + P) - μZ This is a breif summary of the asignment, the whole thing can otherwise be found in proj1mat_eng.pdf and in the plankton.pdf github)

% File for main program (ah_alg)
clear all
close all 
clc

h = 0.5;                  % Step length
h2 = h/2
tspan = 0:h:400;       %Start value, time expressed as days

% Initial conditions
p0 = 20
z0 = 5
x0 = [p0,z0]'; % Write everything in vector form
r = 0.3
K = 108
R = 0.7


%Parameters
alpha = 5.7
mu = 0.024
gamma = 0.05

% Solving Runge knutta 4

X(:,1) = x0;
xin = x0;
xin2 = x0;
Err = [0; 0]; % Matrix for saving error for p and x

time = tspan(1);

for i=1:tspan(end)/h
    
    pout = rk4(@(t,x)lorenz(t,x,K,R,r,alpha,mu,gamma),h,time,xin);
    pout2 = rk4(@(t,x)lorenz(t,x,K,R,r,alpha,mu,gamma),h2,time,xin2);
    X = [X pout];
    if time <= 200
      Err(1,i) = abs(pout2(1)-pout(1));  %Error for p
      Err(2,i) = abs(pout2(2)-pout(2)); %Error for z
    end
    xin = pout;
    xin2 = pout2;
    time = time+h;
end
pout
Err(:,(201:400))= [];
length(Err(1,:))
%length(tspan)
plot(log(tspan(1,(1:201))),log(Err(1,:))) % Erf for p
xlabel('logtime');
%hold on
%plot(log(tspan(1,:)),log(Err(2,:)))
%hold on
%plot(tspan, X(1,:));
%hold on
%plot(tspan,X(2,:))

% Lorenz function
function dp = lorenz(t,x,K,R,r,alpha,mu,gamma)
  dp = [
  r*x(1)*(1-x(1)/K) - R*x(2)*(x(1)^2/(alpha^2 + x(1)^2));
  gamma*R*x(2)*(x(1)^2/(alpha^2 + x(1)^2)) - mu*x(2);
  ];

% Rk4
function pout = rk4(f,h,tk,xk)
%Rungeknutta 4 method
%   This function calculates Rk4
% Startvalue for k1, t = t0

  k1 = f(tk,xk);  
  k2 = f(tk+h/2,xk+h*k1/2);
  k3 = f(tk+h/2,xk+h*k2/2);
  k4 = f(tk+h,xk+h*k3);
  pout = xk + (h/6)*(k1+k2+k3+k4);

I have managed to write the code in ah_alg.m, the function in lorenz and the rk4 in rk4 files. My code produces the correct result for a), but not for b) Can anyone look through my code files and see whats wrong or is it better to start over?