a problem to get a differential equation

Hi, everyone! I have two equations with three variables (lambda_H2O, lambda_H3O and a_H2O), at first I just use a_H2O with 'fsolve' to calculate lambda_H2O and lambda_H3O (so lambda_H2O could be treated as a equation of lambda_H2O(a_H2O).

eq1 = exp(phi_2 * lambda_H2O + phi_1 * lambda_H3O) * lambda_H3O/(1 - lambda_H3O)/(lambda_H2O - lambda_H3O)- K_1 ==0;

eq2 = exp(phi_3 * lambda_H2O + phi_2 * lambda_H3O) * K_2 * (lambda_H2O - lambda_H3O) - a_H2O == 0;

Now I need a differential equation of 'dlambda_H2O / da_H2O', I think lambda_H3O could be represented with lambda_H2O in 'eq1', then substitute it into eq2 and get a relationship of lambda_H2O and a_H2O. I used 'solve', but it failed at first step : (, it shows "Unable to find explicit solution. For options, see help. " Is anyone could tell me how to deal with this problem? Thanks a lot!

The first code is calculation of lambda_H2O with 'fsolve':

function [lambda_H2O, lambda_H3O] = lambda_v(a_H2O,T)
%% calculation of water content in vapor-equilibrium
% Fit parameters for the chemical model of water uptake from the vapor phase
% pre-input values
% a_H2O = preinput.a_H2O; % avtivity of water
% K_1 = preinput.K_1;     % equilibrium coefficient K_1
% T = preinput.T_mem;         % operating temperature in the membrane

E_00v = -41.7;       % binary interaction parameter between H2O and H2O used in the chemical model, g/mol
E_hy0v = -52.0;      % binary interaction parameter between H2O and hydronium ion used in the chemical model, g/mol
E_hyprotv = -3721.6; % binary interaction parameter between hydronium ion and proton used in the chemical model, g/mol
m_mdry = 1100;       % membrane dry equivalent weight (g/mol) [Pukrushpan.Modeling and Control of Fuel Cell Systems and fuel processors] P48 Ta.4.1
R_u = 8.3145;        % universal gas constant (J/mol/K)
K_1 = 100;           % calculation of equilibrium coefficient K_2; -Weber.II
K_2 = 0.217 .* exp(1000./R_u.*(1/303.15-1./T));  % calculation of equilibrium coefficient K_2
phi_1 = 2*(E_00v - 2*E_hyprotv - 2*E_hy0v)/m_mdry;
phi_2 = 2*(E_hy0v - E_00v)/m_mdry;
phi_3 = 2*E_00v/m_mdry;

lambda_ini = [1.2 0.8]; % lambda_ini(1) = lambda_H2O; lambda_ini(2) = lambda_H3O+

lambda = zeros(length (a_H2O),2*length(K_2));
for i = 1 : length (a_H2O)
    % Check if there are multiple variables in K_1 and K_2
    k = length(K_2);
    for j = 1 : k
        [lambda(i, 2*j-1 : 2*j),fval,exitflag] = fsolve(@(x) fun(x),lambda_ini,optimset('Display','off'));
        lambda_H2O(i,j) = lambda(i, 2*j-1); lambda_H3O(i,j) = lambda(i, 2*j);
        ii = 0;
        if exitflag ~= 1
            if ii == 0
                fprintf('non-converged by iteration of lambda! \n')
                ii = ii +1;
            end
            formatSpec = 'error = %e; by a_H2O = %0.3f and T = %3.2f K;\n';
            if length(T) ==1
                fprintf(formatSpec,fval(end),a_H2O(i),T);
            else
                fprintf(formatSpec,fval(end),a_H2O(i),T(j));
            end
        end
    end
end

lambda_H2O = correc_lambda(lambda_H2O);

%% correction of lambda
    function lambda = correc_lambda(x)
        lambda = x.*(1 + exp(0.3 - x));
    end

%% check if K_1 and T are variables
    function F = fun(x)
        if length(K_2) > 1 && length(K_1) == 1
            F = [exp(phi_3 * x(1) + phi_2 * x(2)) * K_2(j) * (x(1) - x(2)) - a_H2O(i);...
                exp(phi_2 * x(1) + phi_1 * x(2)) * x(2)/(1 - x(2))/(x(1) - x(2))- K_1];
        elseif length(a_H2O) >= 1 && length(K_1) == 1 && length(K_2) == 1
            F = [exp(phi_3 * x(1) + phi_2 * x(2)) * K_2 * (x(1) - x(2)) - a_H2O(i);...
                exp(phi_2 * x(1) + phi_1 * x(2)) * x(2)/(1 - x(2))/(x(1) - x(2))- K_1];
        else
            error('issue of setting length in lambda_v.m file!\n');
        end
    end

end

The second is what confusing me now, to get the dlambda_H2O / da_H2O

T = 353.15;
E_00v = -41.7;       % binary interaction parameter between H2O and H2O used in the chemical model, g/mol
E_hy0v = -52.0;      % binary interaction parameter between H2O and hydronium ion used in the chemical model, g/mol
E_hyprotv = -3721.6; % binary interaction parameter between hydronium ion and proton used in the chemical model, g/mol
m_mdry = 1100;       % membrane dry equivalent weight (g/mol) [Pukrushpan.Modeling and Control of Fuel Cell Systems and fuel processors] P48 Ta.4.1
R_u = 8.3145;        % universal gas constant (J/mol/K)
K_1 = 100;           % calculation of equilibrium coefficient K_2; -Weber.II
K_2 = 0.217 .* exp(1000./R_u.*(1/303.15-1./T));  % calculation of equilibrium coefficient K_2
phi_1 = 2*(E_00v - 2*E_hyprotv - 2*E_hy0v)/m_mdry;
phi_2 = 2*(E_hy0v - E_00v)/m_mdry;
phi_3 = 2*E_00v/m_mdry;

syms lambda3 lambda a0

eqn = exp(phi_2 * lambda + phi_1 * lambda3) * lambda3/(1 - lambda3)/(lambda - lambda3)- K_1 ==0;
eqn = rewrite(eqn,'log');
sol = solve(eqn, lambda3,'IgnoreAnalyticConstraints',1)