index exceeds number of array elements (1)

epsilonb = 0.35;
rhoA = 2.057;
lamdaA = 6.5e-1;
Uz = 0.42;
L = 0.16;
Cpa = 1012;
Cpv = 2030;
Y1 = 0.012;
km = 1.5e-3;
a = 0.08;
Daeff = 5e-5;
Yinitial = 0.012;
Tbinitial = 100;
Tps = 25;
Xstar = 0.143;
phigrain = 0.7;
tfinal = 0.010;
deltaz = 0.016;
zinitial = 0;
pW = 21.6236;
Ystar = 0.01264;

n = (L-zinitial)/deltaz + 1;
Y0 = Yinitial*ones(n-1,1);
Tb0 = Tbinitial*ones(n-1,1);
s0 = Y0;
[t,s] = ode23s(@(t,s) gasphasefun(t,s), [0 tfinal],s0);

figure('units','normalized','outerposition',[0 0 1 1]);
plot(t,Y,'-o');

function f = gasphasefun(t,s) 
    epsilonb = 0.35;
    rhoA = 2.057;
    lamdaA = 6.5e-1;
    Uz = 0.42;
    L = 0.16;
    Cpa = 1012;
    Cpv = 2030;
    Y1 = 0.012;
    km = 1.5e-3;
    a = 0.08;
    Daeff = 5e-5;
    Yinitial = 0.012;
    Tbinitial = 100;
    Tps = 25;
    Xstar = 0.143;
    phigrain = 0.7;
    tfinal = 0.010;
    deltaz = 0.016;
    zinitial = 0;
    pW = 21.6236;
    Ystar = 0.01264;
    n = (L-zinitial)/deltaz + 1;
    y = s(2:n-1,1);
    Y = [Yinitial; y];
    dYdt = zeros(n,1);
    dTbdt = zeros(n,1);
    dYdt(1,1) = 0;
    for i = 2:n-1 
        Y = Y(i,1);
        dYdt(i,1) =  Daeff*((Y(i+1)-2*Y(i)+Y(i-1))/deltaz^2)- Uz*((Y(i+1)-Y(i-1))/(2*deltaz)) + a*((1-epsilonb)/epsilonb)*km*(Ystar-Y(i));
    end
    f = dYdt(2:n);
end

[Abolfazl Chaman Motlagh]

please use tool for inserting code or write your code correctly so your question becomes readable.

[DGM]

You convert Y to a scalar

Y = Y(i,1);

and then you try to index into its second and third elements which don't exist:

(Y(i+1)-2*Y(i)+Y(i-1))/deltaz^2

You'll have to decide what do to with Y such that the indexing operations make sense. Perhaps using subscripted operations would simplify this, but I'm not going to dig that deep.