How to extract variable from function file, and how to check is that result real?

I am using matlab function file to solve system of differential equations, and I have variable R in my function file.

function f=fun_z(z,p)
ri=2;
R=ri-z*(ri-1);
sig=1; beta=1;
f=zeros(4,1);
f(1)=-32*1*beta/(R.^4*p(1));
f(2)=(-(2-sig)*8*f(1)/(sig.*R)-f(1)*p(2))/p(1);
f(3)=(-p(2)*f(2)+(2-sig)*(-8*f(2)/R-8*f(1)/(R.*R*p(1)))/sig-f(1)*p(3))/p(1);
f(4)=(-f(2)*p(3)-f(3)*p(2)+(2-sig)*8*(-f(3)/R-f(2)/p(1)+p(2)*f(1)/(p(1).*p(1)))/(sig.*R.*R) -f(1)*p(4))/p(1);

I am not sure did I properly tell to Matlab what is necessary to do with variable R so I extracted it from function file into workspace with command:

[zv,pv,R]=ode45(@fun_z,[1 0],[1; 0; 0; 0])

where I previous in function file changed the first line in:

   function f=fun_z(z,p,R)

and added line of code in function file:

global R;

As a result I need to get R where it is dependent on z, so it means R need to have the same number of elements as z, but that is not the case. I cannot see why, but I only get 6 variables? What would be correct way to tell Matlab exactly what R (R is radial coordinate, z is longitudinal coordinate, R is linearly dependent on z, where z is going from 1 to 0, with some step) is? Or how to extract it properly?

Answer 1

[Jan]

 [zv,pv,R]=ode45(@fun_z,[1 0],[1; 0; 0; 0])

This is a bold and wrong guess. Why do you assume, that ode45 replies an arbitrary global variable from the function to be integrated as 3rd output?

You can either modify the function to accept vectors, then:

function [f, R] = fun_z(z, p)
ri   = 2;
R    = ri - z * (ri - 1);
sig  = 1;
beta = 1;
f       = zeros(4, size(p, 2));
f(1, :) = -32 * beta / (R .^ 4 * p(1, :));
f(2, :) = ((-2+sig) * 8 * f(1, :) / (sig .* R) - f(1, :) *. p(2, :)) / p(1, :);
... etc. I cannot test this currently, please debug this by your own.

Then:

[zv, pv] = ode45(@fun_z, [1 0], [1; 0; 0; 0]);
[~, R]   = fun_z(zv.', pv.');

Or in your case simply:

[zv, pv] = ode45(@fun_z, [1 0], [1; 0; 0; 0]);
ri   = 2;
R    = ri - zv * (ri - 1);