Your ‘melanie_ode_fj’ function has only one output (as it should):
function dx=melanie_ode_fj(t,x)
If you want other outputs from it, create a second version of it (with a slightly different name, I call it ‘melanie_ode_fj_more_outputs’ here) with the outputs you want.
Then, use a for loop and the results of the original ‘melanie_ode_fj’ integration to get the other information you want.
For example (assuming your time vector is ‘t’):
for k = 1: numel(t)
[dx{k} phee{k}] = melanie_ode_fj_more_outputs(t(k),Y(k));
end
NOTE: This is UNTESTED CODE and is for illustration purposes only!
