help with implicit method

Hi guys, This is my first code that I have ever written. I have wrote this code, and I took almost three months to get to this stage. the code is running good but the results are very high, where the output temperature of the material is 250C. this code has a subroutine code where it is consider as an input (Gt). I have changed all the material properties, but nothing change in the output, when I divided Gt (subroutine code)by 10, the output temperature has been divided by 10 and became 25 C, which is acceptable. there is nothing to do with the subroutine, because it is absolutely right.

Is there anybody could help with that please, and find out why my output temperature is very high.

clear all
clc

ho=50;
hi=10;
N=27; %Number of space intervals
M=3600; %number of time steps
Th=196; %thickness of wall in meter;

t_0=0; %initial time
t_f=3600; %final time (sec)
To=34; %exterior wall temperature
k_b=0.7; %thermal conductivity of bricks
roh_b=1600; %density of bricks (kg/m^3)
Cp_b=840; %Cp of bricks (J/kg-K)
T_in=22; %Inner wall temperature
abs=0.7;
%dx = (Th/(N-1))*0.001; %size of the mesh
 dx=0.002;
x = 0:dx:N-1;
x = x'; %inverse of mesh
bi=hi*dx/k_b; %coefficient of terms
bo=ho*dx/k_b;

%n-Octadecane PCM
k_ps=0.358;
k_pl=0.148;
Cp_ps=1.934;
Cp_pl=2.196;
roh_ps=865;
roh_pl=780;
lamda=2.435;
dT=0.2;

Tm=5.5;

Gt=solarradiationcode;
dz=2;

xr = length(Gt);
Te = zeros(xr,1);
Final(1,xr) = zeros;

dt =3;
%initial condition
u(:,1)=ones(N,1)*15;

    for h = 1:length(Gt)
        Te(h) = To*(abs*Gt(h)/ho);
      for j=1:1:M  %time step

    alfa= k_b/(Cp_b*roh_b);
      % define the ratio r

      r=alfa*dt/dx^2; %Fourier number

      Beta=[2*r*bo*Te(h); zeros(N-2,1); 2*r*bi*T_in];

      A = zeros(N,N);
      A(1,1) = 1+2*r+2*r*bo;
      A(1,2) = -2*r;

       for i=2:1:N-1
  %         
  %         
          if i<N/3 || i>2*N/3
             '1'
              alfa= k_b/(Cp_b*roh_b);
              r=alfa*dt/dx^2;
          end
  %         
          if i==N/3 || i==2*N/3
  % %        '2'
  %             
          Cp=2*Cp_b*Cp_ps/(Cp_ps+Cp_b);

              roh=2*roh_ps*roh_b/(roh_ps+roh_b);

              k=2*k_b*k_ps/(k_ps+k_b);

              alfa= k/(Cp*roh);
              r=alfa*dt/dx^2;
          end

          if i>N/3 && i<2*N/3
  %            '3'
              T=u(i,j);
              if T<Tm

                  Cp=Cp_ps;
                  roh=roh_ps;
                  k=k_ps;

               elseif T>Tm-dT && T<Tm+dT

                  Cp= ((Cp_ps+Cp_pl)/2)+((roh_ps+roh_pl)/2)*(lamda/dT);
                  roh=roh_ps;
                  k=k_ps;

              elseif T>Tm+dT

                  Cp=Cp_pl;
                  roh=roh_pl;
                  k=k_pl;

              end             
              alfa=k/(Cp*roh);
              r=alfa*dt/dx^2;
          end

          A(i,i-1) = -r;
          A(i,i) = 1+2*r;
          A(i,i+1) = -r;
      end

      alfa= k_b/(Cp_b*roh_b);
      % define the ratio r
      r=alfa*dt/dx^2; %Fourier number

      %A(N,N+1) = 0;
      A(N,N-1) = -2*r;
      A(N,N) = 1+2*r+2*r*bi;

      [L,UU] = LU(A);

      % note: matlab can tell the difference between u and U, so I could have U as the
      % matrix and u as the solution.  But not all languages will do this, sometimes you
      % have to tell them explicitly to care about case.  The safest thing is to not have the
      % habit of depending on cases.

      % then we'll just do
      %[y] = down_solve(L,u_old);
      %[u_new] = up_solve(UU,y);
      b=u(:,j)+Beta;

      [y] = down_solve(L,b);
      u(:,j+1) = up_solve(UU,y);

      end
      Final(1,h) = u(15,M);
      % Final(1,h) = u(15,M);
      %Final(1,h) = u(15,M);

    end
  plot(Final)

[George Papazafeiropoulos]

What are the functions down_solve and up_solve?

[Libya]

thanks George for your reply.

I have an A matrix at each time step, So I have used a matrix elimination to solve the A matrix. The LU(A), down_solve and up_solve are a sub-code for matrix elimination method.

[Image Analyst]

Well, attach them. Attach all needed m-files with the paper clip icon so we can try your code. Also, read this: http://www.mathworks.com/matlabcentral/answers/13205-tutorial-how-to-format-your-question-with-markup

[Libya]

When I went through this code line by line I found that, when I eliminate the A matrices at each time step the output increases. Anybody could help.