I am trying to solve ODE having more than 1 dependent variable. but I am not able to solve it with dsolve function. please help me to find correct function to solve ODE.This is an mechanical engineering equation so it is bit complicated and lengthy

syms Temp_w Temp_g Temp_col
%defining constant
m_w = 15; % mass of water in Kg
m_g = 5; % mass of glass cover in Kg
m_col = 6; % mass of collecter in Kg

c_w = 4179.6; % specific heat of water
c_g = 800; % specific heat of glass
c_col = 2000; % specific heat of collecter

a_w = 0.05; % absorbtivity of water
a_g = 0.08; % absorbtivity of glass
a_col = 0.9; % absorbtivity of collecter

Temp_a = 300; % ambient temerature in Kelvin

G_max = 630; % maximum irradiation of sun
G = 1:10;
for time = 1:10
        G(time) = (G_max/2)*(sin((pi*time)/12));
end
G_avg = mean(G); % avg of solar irradiation during a day

% belove are term use in ODE equation
q_ev = 16.273 * 0.884 * (((Temp_w - Temp_g + (((11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13))))*(Temp_w))/(268900 - 11.6834 - (3816.44/(Temp_w-46.13)))))^(1/3))) * (11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13)))) * (11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13))));

q_cw = 0.884 * (((Temp_w - Temp_g + (((11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13))))*(Temp_w))/(268900 - 11.6834 - (3816.44/(Temp_w-46.13)))))^(1/3))) * (11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13)))) * (Temp_w - Temp_g);

q_rw = 5.67 * (10^(-8)) * 0.96 * ((Temp_w^4) - (Temp_g^4));

q_ca = 6.4 * (Temp_g - Temp_a);

q_ra = 5.67 * (10^(-8)) * 0.93 * ((Temp_g^4) - (Temp_a^4));

q_w = 1000 * (Temp_col - Temp_w);

q_ins = (0.026 * (Temp_col - Temp_a))/0.01;

ag_G = a_g*G_avg;
aw_G = a_w*G_avg;
acol_G = a_col*G_avg;

% ODE equation(1) = q_ev + q_cw + q_rw + ag_G - q_ra - q_ca == m_g*c_g* (dTemp_g/dt)

syms Temp_g(t)
ode_1 = m_g*c_g*diff(Temp_g,t) == 16.273 * 0.884 * (((Temp_w - Temp_g + (((11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13))))*(Temp_w))/(268900 - 11.6834 - (3816.44/(Temp_w-46.13)))))^(1/3))) * (11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13)))) * (11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13)))) + 0.884 * (((Temp_w - Temp_g + (((11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13))))*(Temp_w))/(268900 - 11.6834 - (3816.44/(Temp_w-46.13)))))^(1/3))) * (11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13)))) * (Temp_w - Temp_g) + 5.67 * (10^(-8)) * 0.96 * ((Temp_w^4) - (Temp_g^4)) + (a_g*G_avg) - 6.4 * (Temp_g - Temp_a) - 5.67 * (10^(-8)) * 0.93 * ((Temp_g^4) - (Temp_a^4));
Temp_g_sol(t) = dsolve(ode_1)

% ODE equation(2) = -q_ev - q_cw - q_rw + aw_G + q_w == m_w*c_w* (dTemp_w/dt)

syms Temp_w(t)
ode_2 = m_w*c_w*diff(Temp_w,t) == -(16.273 * 0.884 * (((Temp_w - Temp_g + (((11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13))))*(Temp_w))/(268900 - 11.6834 - (3816.44/(Temp_w-46.13)))))^(1/3))) * (11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13)))) * (11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13))))) - (0.884 * (((Temp_w - Temp_g + (((11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13))))*(Temp_w))/(268900 - 11.6834 - (3816.44/(Temp_w-46.13)))))^(1/3))) * (11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13)))) * (Temp_w - Temp_g)) - (5.67 * (10^(-8)) * 0.96 * ((Temp_w^4) - (Temp_g^4))) + (a_w*G_avg) + (1000 * (Temp_col - Temp_w));
Temp_w_sol(t) = dsolve(ode_2)

% ODE equation(3) = acol_G - q_w - q_ins == m_col*c_col* (dTemp_col/dt)

syms Temp_col(t)
ode_3 = m_col*c_col*diff(Temp_col,t) == (a_col*G_avg) - (1000 * (Temp_col - Temp_w)) - ((0.026 * (Temp_col - Temp_a))/0.01);
Temp_col_sol(t) = dsolve(ode_3)

[Alan Stevens]

It looks like you only have one independent variable (time). Temp_w is a dependent variable. Do you know the values of Temp_g and Temp_a? If so, specify them. If they are further dependent vaiables. then you need more equations (differential or algebraic) to allow them to be determined (simultaneously with Temp_w if necessary).

[RUTVIK RAVAL]

@Alan Stevens, you are right. actully there are 3 equation and 3 unknowns. But I thought if I know how to solve 1 equation then I can solve other 2 also. But as you mentioned i have edited the code so that you might tell me the way to solve it. Actully my aim to solve all 3 ode so I get 3 equation and after solving that 3 equation I may get the value of Temp_w, Temp_g and Temp_col. Please help me in code if you know how to solve above 3 ode and then that 3 equations. Thank you very much @Alan Stevens

[RUTVIK RAVAL]

@Alan Stevens. in solving ODEs if you need to give initial condition than give any constant value

Answer 1

[Alan Stevens]

Here's a possible way using ode45. I've used arbitrary initial values for the three temperatures, so you will need to replace them with the true values.

a_w = 0.05; % absorbtivity of water
a_g = 0.08; % absorbtivity of glass
a_col = 0.9; % absorbtivity of collecter

G_max = 630; % maximum irradiation of sun
G = 1:10;
for time = 1:10
        G(time) = (G_max/2)*(sin((pi*time)/12));
end
G_avg = mean(G); % avg of solar irradiation during a day

ag_G = a_g*G_avg;
aw_G = a_w*G_avg;
acol_G = a_col*G_avg;
absorp = [ag_G, aw_G, acol_G];

tspan = [1 10];
Temp_g0 = 310;    %%% ?   Replace with true values.
Temp_w0 = 320;    %%% ?
Temp_col0 = 330;  %%% ?
% Combine the three initial temperatures into a single vector.
T0 = [Temp_g0, Temp_w0, Temp_col0];

% Call the function that calculates the rates of change with ode45
% using the following syntax (try doc ode45 for more details).
[t, T] = ode45(@(t,T) fn(t,T,absorp),tspan,T0);

% Extract the individual values of the three temperatures from T.
Temp_g = T(:,1); Temp_w = T(:,2); Temp_col = T(:,3);

% Display the results
plot(t,Temp_g,t,Temp_w,t,Temp_col),grid
xlabel('time'),ylabel('Temps')
legend('Tg','Tw','Tcol')

% The rate of change of temperatures wrt time are contained in the following function. 
% The input arguments must be t and T or ~ and T (the latter if t is not used explicitly 
% within the function). Also, any extra data that is used within the function
% ca be passed in as an extra argument - like absorp here.  
% The three temperatures are passed in via a single vector T.
% Their individual rates of change are returned as a column vector in dTdt. 
function dTdt = fn(~,T,absorp)
        
        % Extract the three individual temperatures for use in the rates 
        % of change equations below.
        Temp_g = T(1); Temp_w = T(2); Temp_col = T(3);
        
        % Extract the individual pieces of data passed to the function
        % in absorp.
        ag_G = absorp(1);
        aw_G = absorp(2);
        acol_G = absorp(3);
        
        %defining constant
        m_w = 15; % mass of water in Kg
        m_g = 5; % mass of glass cover in Kg
        m_col = 6; % mass of collecter in Kg

        c_w = 4179.6; % specific heat of water
        c_g = 800; % specific heat of glass
        c_col = 2000; % specific heat of collecter

        Temp_a = 300; % ambient temerature in Kelvin

        % ODE equation(1) = q_ev + q_cw + q_rw + ag_G - q_ra - q_ca == m_g*c_g* (dTemp_g/dt)
        dTemp_gdt = (16.273 * 0.884 * (((Temp_w - Temp_g + ...
            (((11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - ...
            (3816.44/(Temp_g-46.13))))*(Temp_w))/(268900 - 11.6834 - ...
            (3816.44/(Temp_w-46.13)))))^(1/3))) * (11.6834 - ...
            (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13)))) *...
            (11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13)))) +...
            0.884 * (((Temp_w - Temp_g + (((11.6834 - (3816.44/(Temp_w-46.13)) - ...
            (11.6834 - (3816.44/(Temp_g-46.13))))*(Temp_w))/(268900 - 11.6834 - ...
            (3816.44/(Temp_w-46.13)))))^(1/3))) * (11.6834 - (3816.44/(Temp_w-46.13)) - ...
            (11.6834 - (3816.44/(Temp_g-46.13)))) * (Temp_w - Temp_g) + 5.67 * (10^(-8)) * ...
            0.96 * ((Temp_w^4) - (Temp_g^4)) + ag_G - 6.4 * (Temp_g - Temp_a) - ...
            5.67 * (10^(-8)) * 0.93 * ((Temp_g^4) - (Temp_a^4)))/(m_g*c_g);

        % ODE equation(2) = -q_ev - q_cw - q_rw + aw_G + q_w == m_w*c_w* (dTemp_w/dt)
        dTemp_wdt = (-(16.273 * 0.884 * (((Temp_w - Temp_g + (((11.6834 - ...
            (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13))))*...
            (Temp_w))/(268900 - 11.6834 - (3816.44/(Temp_w-46.13)))))^(1/3))) *...
            (11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13)))) *...
            (11.6834 - (3816.44/(Temp_w-46.13)) - (11.6834 - (3816.44/(Temp_g-46.13))))) - ...
            (0.884 * (((Temp_w - Temp_g + (((11.6834 - (3816.44/(Temp_w-46.13)) - ...
            (11.6834 - (3816.44/(Temp_g-46.13))))*(Temp_w))/(268900 - 11.6834 - ...
            (3816.44/(Temp_w-46.13)))))^(1/3))) * (11.6834 - (3816.44/(Temp_w-46.13)) - ...
            (11.6834 - (3816.44/(Temp_g-46.13)))) * (Temp_w - Temp_g)) - (5.67 * (10^(-8)) * ...
            0.96 * ((Temp_w^4) - (Temp_g^4))) + aw_G + (1000 * (Temp_col - Temp_w)))/(m_w*c_w);

        % ODE equation(3) = acol_G - q_w - q_ins == m_col*c_col* (dTemp_col/dt)
        dTemp_coldt = (acol_G - (1000 * (Temp_col - Temp_w)) - ...
            ((0.026 * (Temp_col - Temp_a))/0.01))/(m_col*c_col);
       
        % Assemble a column vector of the values of the rates of change.
        dTdt = [dTemp_gdt; dTemp_wdt; dTemp_coldt]; 

end

[RUTVIK RAVAL]

Thank you very much @Alan Stevens for your efforts. I am currently doing my last year project in final year in mechanical engineering. I was not going further in my project because of meshed up this code. But you help me and it will very useful to me. Again thank you very much for quick reply and solving my matlab problem. Still I request you to put some comments in the code which will help me to understand this code better. I really appritiate your work

[Alan Stevens]

Ok. I've added a few comments. You should also study the documentation on ode45.