Extract value from function

Hello friends, I have this code that return values for dc1dt and dq1dt. But I want to see the values for dc1dz, how can I do that?

function dydt=f(t,y,R,P,epsb,u,LDF1,b1,qm1,c10,Nz,dz,rhop,nco2,T)
dydt=zeros(length(y),1);
dc1dt=zeros(Nz,1);
dq1dt=zeros(Nz,1);

%Assign values
c1=y(1:Nz);
q1=y(Nz+1:2*Nz);

%BC
c1(1)=c10;
c1(end)=(4*c1(end-1)-c1(end-2))./3;

%interior
for i=2:Nz-1
    
    dc1dz(i)=(c1(i+1)-c1(i-1))./(2.*dz);
    q1star(i)=qm1.*((b1.*0.197).^(1/nco2))./(1+(b1.*0.197).^(1/nco2)); 
    dq1dt(i)=LDF1.*(q1star(i)-q1(i));
    
   
    %main f
    dc1dt(i)=-u*dc1dz(i)-(rhop.*R.*T./P *((1-epsb)./epsb).*dq1dt(i));
    

end
dydt=[dc1dt;dq1dt];
end

Answer 1

[Walter Roberson]

Chances are that you do not want to see that.

You are obviously calling an ode*() function. Your function f is going to be evaluated at a series of trial points "near" the current location, in order to predict what the new location should be, and then a cross-check is done to see if the value at another location is "close enough" to a value that the proposed coefficients would predict. For ode45() this involves calculating the function 6 times at each location, with none of the computed locations necessarily making it into the output. If the cross-check was too far away from prediction, then the ode*() goes back and tries again with a smaller step size, another 6 calculations. For any system that cannot be described by a quintic or lower degree polynomial, it is normal for a number of points to be rejected.

Knowing the dc1dz values for a number of trial locations that were never going to make it to output is seldom of benefit.

None of the locations evaluated at necessarily make it to output: the output can end up consisting entirely of predicted values.

For these reasons, most of the time it is much better to call ode45() or whatever normally, getting outputs for various time steps, and then go back and compute the dc1dz values that would be associated with those boundary conditions.

You can use the same function with dc1dz as the second output; loop over the rows of boundary conditions calling the function and capturing the second output... and do whatever with it.