error in ode45 - must return a column vector

Hi everyone,

Im trying to solve a PDE using ode45. But the following error are ocurring: 'EDP must return a column vector'.

I put my code blow for you guys help me to find the error.

------------------------------------------------------------------

global pipe Tin I Tg W T TW;

T = [];

Ts = [];

TW = [];

Tws = [];

% condições iniciais

Iini = IC151(1,2);

Tgini = TA075(1,2);

Wini = FA061(1,2);

Tini = TA072(1,2);

T = Tini;   

pipe.Tref = 250;  

pipe.rho = 903 - 0.672 * pipe.Tref;     

pipe.Cp = 1820 + 3.76 * pipe.Tref;   

pipe.rhow = 7800;      

pipe.Cw = 450;  

pipe.ri = 0.009; 

pipe.ro = 0.013; 

pipe.di = pipe.ri * 2;   

pipe.do = pipe.ro * 2;  

pipe.Ai = pi * pipe.ri^2; 

pipe.Ao = pi * pipe.ro^2; 

pipe.ho = 11;

pipe.hi=800;   

pipe.L=172;   

pipe.S=150; 
                                                       pipe.dX=pipe.L /pipe.S;    

pipe.n0=0.43;  

pipe.G=1.82;  
                                                     pipe.vs=Wini /pipe.Ai;   
pipe.gamma = (pipe.n0 * pipe.G) / (pipe.Ao * pipe.rhow * pipe.Cw);

pipe.Tau1 = (pipe.Ai * pipe.rho * pipe.Cp) / (pi * pipe.di * pipe.hi);

pipe.Tau12 = (pipe.Ao * pipe.rhow * pipe.Cw) / (pi * pipe.di * pipe.hi);

pipe.Tau2 = (pipe.Ao * pipe.rhow * pipe.Cw) / (pi * pipe.do * pipe.ho);

pipe.c = (pipe.vs * pipe.Tau1 * (1 + (pipe.Tau2 / pipe.Tau12)));

TW = (Tini * pipe.Tau2) / (pipe.Tau12 + pipe.Tau2) + (Tgini * pipe.Tau12) / (pipe.Tau12 + pipe.Tau2) + (Iini * pipe.gamma * pipe.Tau12 * pipe.Tau2) / (pipe.Tau12 + pipe.Tau2);

for(i=1:pipe.S)

    T(i+1) = pipe.gamma * pipe.Tau2 * Iini + Tgini + (T(1) - pipe.gamma * pipe.Tau2 * Iini - Tgini) * exp(-(i * pipe.dX) / pipe.c);

    TW(i+1) = (T(i+1) * pipe.Tau2) / (pipe.Tau12 + pipe.Tau2) + (Tgini * pipe.Tau12) / (pipe.Tau12 + pipe.Tau2) + (Iini * pipe.gamma * pipe.Tau12 * pipe.Tau2) / (pipe.Tau12 + pipe.Tau2);

end

Ts = T(1,pipe.S+1);

Tws = TW(1,pipe.S+1);

n = 2;

for k = 1:2481

      I = IC151(n,2);              % from workspace

      Tin = TA072(n,2);            % from workspace

      Tg = TA075(n,2);             % from workspace

      W = FA061(n,2);              % from workspace

      [dT dTW] = ode45(@EDP,[0 1],[T,TW],[]);

      n= n+1; 

end

------------------------------------------------------------------

and the EDP file has:

-------------------------------------------------------------------

function [dT, dTW] = EDP(t,x)

global pipe W I Tg Tin T TW;

T = x(1:pipe.S+1);

T(1) = Tin;

TW = x(pipe.S+1+1:2*(pipe.S+1));

vs = W / pipe.Ai;

dT(1) = 0;

dTW(1) = 0;

for(i=2:pipe.S+1)

    dT(i) = -vs * ((T(i) - T(i-1)) / pipe.dX) + ((TW(i) - T(i)) / pipe.Tau1);

    dTW(i) = pipe.gamma * I + ((Tg - TW(i)) / pipe.Tau2) - ((TW(i) - T(i)) / pipe.Tau12);

end

------------------------------------------------------------------

I would like to you help me.

best regards,

Gustavo

[Jan]

Please format your code to improve the readability. See the "Markup help" link.

Answer 1

[Jan]

Simply replace the "dT(1)=0" line to pre-allocate dT in the correct shape:

dT = zeros(pipe.S+1, 1);

The same for dTW.

[Gustavo]

thanks Jan Simon, now my simulation is running :)

best regards,

Gustavo