solving an equation system with time and space derivation

Hi everyone,

i am trying to solve an equation system with time and space derivation.i used ode45 to solve this system but there is an error messege .could you help me to run the simulation

here is the function code

function  dy= solution(t,z)
global a2 a3 a4 alphaverre rhopv Cppv deltapv APTV  n N   b2 b3 b4 c1 c2 c3 c4 d1 d2 i j  hi dt 
% Set aside storage
Tglass = zeros(N,n); 
Tpvt   = zeros(N,n);
Twax   = zeros(N,n);
Tfluid = zeros(N,n);
%t = zeros(1,n);
%z = zeros(1,N);
% Initialise temperatures
Tglass(1,1)  = 29;
Tpvt(1,1)    = 25;
Twax(1,1)    = 20;
Tfluid(1,1) = 19;
for i = 2:N   % space loop
   Tglass(i,1) = 29;
   Tpvt(i,1)   = 25;
   Twax(i,1)   = 20;
   Tfluid(i,1) = 19;
   z(i) = (i-1).*hi;
   
   for j = 2:n  % time loop
       
       % Use a simple Euler integration for simplicity
       t(j) =(j-1).*dt;
       sinterm = 0.08.*sin(t(j)./6000);
      %sinterm= 0.08*(sin(t(j)/tau-pi/2)+1)*950/0.16+50;
       dTglassdt = 0.105+7.8.*10^-6.*sinterm -a2.*Tglass(i,j-1) -a3.*Tglass(i,j-1)^4 + a4.*Tpvt(i,j-1);
       dTpvtdt = (1-alphaverre).*sinterm./(rhopv*Cppv*deltapv*APTV)- b2.*Tglass(i,j-1) + b3.*Tpvt(i,j-1) + b4.*Twax(i,j-1);
       dTwaxdt =  c1 + c2.*Tpvt(i,j-1) - c3.*Twax(i,j-1) - c4.*Tfluid(i,j-1);
       dTfluiddz = (Tfluid(i,j-1)- Tfluid(i-1,j-1))./hi;
       dTfluiddt = d1.*(Twax(i,j-1) - Tfluid(i,j-1)) - d2.*dTfluiddz;
       
       Tglass(i,j)= Tglass(i,j-1) + dt.*dTglassdt;
       Tpvt(i,j) = Tpvt(i,j-1) + dt.*dTpvtdt;
       Twax(i,j)= Twax(i,j-1) + dt.*dTwaxdt;
       Tfluid(i,j)= Tfluid(i,j-1) + dt.*dTfluiddt;
    end
    
end
dy=[sinterm;Tglass(i,j);Tpvt(i,j);Twax(i,j);Tfluid(i,j)];
end 

and for the script i write the following lines

clc
clear
global a2 a3 a4  b2 b3 b4 c1 c2 c3 c4 d1 d2 N n L tf alphaverre Uv  Ta sigma epsilonv Averre   rho Cp  kverre rhopv Cppv  deltav kPVT APVT deltapv m Cpf Tfin rhocire Cpcire h mf  cpf rhof Acoil d hf AWAX ewax hi dt
alphaverre=0.03;
Uv=0.15; 
Averre=0.8571;
Ta=25;
sigma=5.670367*10^-8;
epsilonv=0.88;
rho=3000;
Cp= 500;
kverre=1;
rhopv=2330;
Cppv=702.9;
deltav=0.003;
kPVT=149;
APVT =0.8571;
deltapv=0.0002;
m=0.12;
Cpf=4180;
Tfin=25;
rhocire=930;
Cpcire=2100;
h=30;
mf=0.12;
cpf=4180;
rhof=1000;
Acoil=6.3*10^-5;
d=2*sqrt(Acoil/pi);
hf=166.216;
AWAX=0.8571;
ewax=0.002;
a2=(((5.7+3.8.*Uv).*Averre)/ rho*Cp*Averre*deltav)+(kverre.* Averre)/( rho*Cp*Averre*deltav*deltav);
a3=(sigma*epsilonv* Averre)/( rho*Cp*Averre*deltav);
a4=(kverre* Averre)/(rho*Cp*Averre*deltav*deltav);
b2=(Averre*kverre)/((rhopv*Cppv*deltapv*APVT)*deltav);
b4=kPVT*APVT/(rhopv*Cppv*deltapv*APVT*deltapv);
c1=((h*APVT*Ta)+(m*Cpf* Tfin))/(rhocire*Cpcire*AWAX*ewax);
c2=(kPVT*APVT)/(rhocire*Cpcire*AWAX *ewax*deltapv);
c3=(kPVT*APVT)/((rhocire*Cpcire*AWAX *ewax)*deltapv)+h*APVT/(rhocire*Cpcire*AWAX *ewax);
c4=(m*Cpf)/(rhocire*Cpcire*AWAX *ewax);
d1=pi*d*hf/rhof*cpf*Acoil;
d2=(mf*cpf)/(cpf*rhof*Acoil);
b3= (kverre*Averre/((rhopv*Cppv*APVT*deltapv)*deltav))-(kPVT*APVT/((rhopv*Cppv*deltapv*APVT)*deltapv));
% Spatial geometry
L =  2 ;      % Insert length 
N =  20 ;      % insert number of temperature nodes
hi = L/(N-1); % length of section (ie distance between temperature nodes
% Timespan
tf = 36000;        % insert final time
n  = 10;       % insert number of timesteps for which you want output
dt = tf/n;     % timestep

[t,z]=ode45(@solution,[1 36000],[70;29;25;20;19]);
figure(1)
plot (t,sinterm,'x')
title('FLUX SOLAIRE');
xlabel('Time t');
ylabel('puissance solaire');
legend('G')
hold on
figure(2)
plot(t,Tglass,'r')
hold on
plot(t,Tpvt,'y')
hold on
plot(t,Twax,'g')
hold on
title('temperature variation in time');
xlabel('Time t');
ylabel('temperature T  °C');
legend('T_g','T_pvt','T_wax','T_fluide')
figure(3)
subplot(2,1,1)
plot(t,Tfluid),grid
xlabel('time'),ylabel('Tfluid')
subplot(2,1,2)
plot(z,Tfluid),grid
xlabel('space'),ylabel('Tfluid')