function [t,x] = rk4m(f,tspan,x0,n),
global conta
conta=0;
a = tspan(1);
b = tspan(2);
h = (b-a)/n;
t = (a+h:h:b);
k1 = h*feval(f,a,x0);
k2 = h*feval(f,a+h/2,x0+k1/2);
k3 = h*feval(f,a+h/2,x0+k2/2);
conta=4;
k4 = h*feval(f,a+h,x0+k3);
conta=0;
x(1,1:2) = x0 + (k1)/6+ (k2)/3 + (k3)/3 + (k4)/6;
for i=1:n-1
k1 = h*feval(f,t(i),x(i,1:2));
k2 = h*feval(f,t(i)+h/2,x(i,1:2)+k1/2);
k3 = h*feval(f,t(i)+h/2,x(i,1:2)+k2/2);
conta =4; %añadido el conta para que solo una vez de las llamadas a la funcion calcule el caudal de salida, la idea
%es independizar la implementacion totalmente del modelo
k4 = h*feval(f,t(i)+h,x(i,1:2)+k3);
conta=0;
x(i+1,1:2) = x(i,1:2) + (k1)/6+ (k2)/3 + (k3)/3 + (k4)/6;
end
conta =4;
k4dummy = h*feval(f,t(n)+h,x(n,1:2)+k3);
conta=0;
t = [a t];
v1 = [x0(1) x(:,1)'];
v2 = [x0(2) x(:,2)'];
x(1:n+1,1)=v1';
x(1:n+1,2)=v2';
