MATLAB: How to reduce the time

vpasolve inside an ode

I'm writing up a code to solve a system of mass balances. The balances are ordinary differential equations, so I used ode23s to solve it. Additionally, inside the function that has the differential equations, I used vpasolve because I have another system of equations for the parameters of one of the differential equation. I'm solving the ode within a loop like this:
 for i=1:5 % Hours
        Mi(i,1)=O(i); Mi(i,2)=C(i); Mi(i,3)=W(i);
        auxt(1)=0; 
        ti=[auxt(i) i];
        [t,M]=ode23s(@TrayFunction,ti,Mi(i,:),[],par);
        auxt(i+1)=auxt(i)+1;
        O(i+1)=M(end,1); C(i+1)=M(end,2); W(i+1)=M(end,3);
        aw(i+1)=aw_v; H(i+1)=H_v; phi(i+1)=phi_v;
    end
And the function is:
    function dM = TrayFunction(t, M, par)
    % Components
        O=M(1);       % Oxygen concentration $\frac{molO_2}{m^3_{a}}$
        C=M(2);       % Carbon dioxide concentration $\frac{molCO_2}{m^3_{a}}$
        W=M(3);       % Total amount of water in a control volume $\frac{kg_{H_2O}}{m^3_{b}}$
    % Complementary equations
        syms phi_v
        % Moisture content in the solids
        eqn2=(1-par.epsilon)*par.rho_s*phi_v+par.epsilon*par.rho_a*((par.Psat/(0.0288*par.P))*(1-exp(-0.0004373*(par.T+57.926)*phi_v^(1.2626))))/(55.5556-((55.5556*par.Psat)/par.P)*(1-exp(-0.0004373*(par.T+57.926)*phi_v^(1.2626))))-W==0;
        phi_v=vpasolve(eqn2, phi_v, par.phi); phi_v=double(phi_v);
        aw_v=isotherm(phi_v,par);
        H_v=definition(aw_v,par);
        assignin('base', 'aw_v', aw_v); assignin('base', 'H_v', H_v); assignin('base', 'phi_v', phi_v)
    % Oxygen balance
        dO=(1/(par.epsilon*par.At*par.dz))*(par.Ky_O*par.epsilon*par.At*(par.O_h-O)-(1-par.epsilon)*par.rho_s*par.At*par.dz*par.r_O);
    % Carbon dioxide balance
        dC=(1/(par.epsilon*par.At*par.dz))*((1-par.epsilon)*par.rho_s*par.At*par.dz*par.r_C-par.Ky_C*par.epsilon*par.At*(C-par.C_h));
    % Water balance 
        dW=(1/(par.At*par.dz))*(-par.Ky_eW*par.At*(aw_v-par.aw_eq)-par.Ky_dW*par.epsilon*par.At*par.rho_a*(H_v-par.H_h));
        dM=[dO dC dW]';
    end
Does anyone have any suggestion in order to have a more efficient code? It is taking too much time. I've thought that the problem was in the use of assignin, but it is too slow even if I remove that part.

Best Answer

Get rid of the sym. Change eqn2 into a function handle with phi_v as parameter. Change the vpasolve into an fsolve(eqn2, par.phi)
Related Question