How do I fix this runge-kutta algorithm implimentation?

I have a matrix c with 3 columns containing 3 sets of data I want to plot, mimicing 4 bolus injections total, changing the "q" value to modify blood concentration at each injection to maintain a certain concentration level. I used the following documentation as a guide to establish the algorithm itself but I have reached a standstill: https://www.mathworks.com/videos/solving-odes-in-matlab-3-classical-runge-kutta-ode4-117528.html

%4 bolues total - at time 0, 24, 48, 72
%1st function - dcdt=f-mc
%this is dydt in the documentation
dcdt=@(t,c,f)(f-Mc);
%impliment algorithm
%pick parameters of start and stop time and step size
%solve for what f should be in each step of time
%use create_f, give it t, h, q, feed it into dcdt to calculate k values
%update c and t
t0=0; %start time
tfinal=94; %stopped time
h=0.1; %step size
%initial conditions
c=1;
c0=0;
t=(0:h:tfinal)';
%dcdt=create_f(t,h,5)
something=ode4(dcdt,t0,h,tfinal,c0);
plot(t,something)
%runge-kutta algorithm implimentation
function cout=ode4(dcdt,t0,h,tfinal,c0)
c=c0;
cout=c;
for t=t0:h:tfinal-h 
    k1=dcdt(t, c);
    k2=dcdt(t+(h/2), c+h*(k1/2));
    k3=dcdt(t+(h/2), c+h*(k2/2));
    k4=dcdt(t+h, c+h*k3);
    c=c+h*(k1+2*k2+2*k3+k4)/6;
    cout=[cout; c];
end
end
%2nd function - define function to create f based on time
%f = [f1; f2; f3], f2=f3=0 for all time
%f1=sum of boluses
%must take time, step size, q which is delta(t-tb), q is plasma []
%immediately after bolus injuction. play with q until i get the target concentrations within the bounds
function out=create_f(t,h,q)
f1=(q/h)*((u(t)-u(t-h))+(u(t-24)-u(t-24-h))+(u(t-48)-u(t-48-h))+(u(t-72)-u(t-72-h)));
f2=0;
f3=0;
out=[f1; f2; f3];
end
%3rd function - unit step, takes time, returns output of unit step function out=u(t)
if t>=0
    out=1;
else 
    out=0;
end
end

[Jan]

Please explain, what "I have reached a standstill" means. You do know, what the problem is, but the readers don't. Instead of spending time to help you, they have to guess at first, what you consider as problem.

[Torsten]

dcdt=@(t,c,f)(f-Mc);

M is not defined.

And most probably you mean M*c instead of Mc.

k1=dcdt(t, c);
k2=dcdt(t+(h/2), c+h*(k1/2));
k3=dcdt(t+(h/2), c+h*(k2/2));
k4=dcdt(t+h, c+h*k3);

It will work, but you should notice that dcdt is defined with three arguments (t,c,f).