I don't see any reason why I should be exceeding matrix dimensions here. The constants in the function (z onwards till Dinf) are defined in a separate script,which I 'run' just before executing the function.

`function [solution]=nano2(x,z,PHI,gamma,Kc,Dp,PHI_B,Mc,N,N_c,deltax_j,J_v,f,gc,T,Dinf)A=x(1:3,1:N+2); %Number of rows to be modified with number of components`

B=x(4:6,1:N+2);Q=x(7:9,1:N+2);R=x(10:12,1:N+2);E=x(13:15,1:N+2);F=x(16:18,1:N+3); %Might need to look at its dimensions.

PSI=x(19,1:N+3);c=x(20:22,1:N+3);zeta=x(23,1);for c1=1:N_cfor c2=3:N+2solution=[A(c1,c2)-((Dp(c1)/deltax_j)-(1/2)*Kc(c1)*J_v+(1/2)*z(c1)*Dp(c1)*(f/(gc*T))*((PSI(c2+1)-PSI(c2))/deltax_j));%Membrane active layer domain for ion i-linearization of Nernst-Planck equation--------(25-31)

B(c1,c2)+(Dp(c1)/deltax_j)-(1/2)*Kc(c1)*J_v+(1/2)*z(c1)*Dp(c1)*(f/(gc*T))*((PSI(c2+1)-PSI(c2))/deltax_j);Q(c1,c2)-(1/2)*z(c1)*Dp(c1)*(f/(gc*T))*((c(c1,c2+1)+c(c1,c2))/deltax_j);R(c1,c2)+(1/2)*z(c1)*Dp(c1)*((f/(gc*T)))*((c(c1,c2+1)+c(c1,c2))/deltax_j);E(c1,c2)+J_v;F(c1,c2)+(1/2)*z(c1)*Dp(c1)*(f/(gc*T))*(c(c1,c2)+c(c1,c2+1))*((PSI(c2+1)-PSI(c2))/deltax_j);F(c1,c2)-A(c1,c2)*c(c1,c2)-B(c1,c2)*c(c1,c2+1)-Q(c1,c2)*PSI(c2)-R(c1,c2)*PSI(c2+1)-E(c1,c2)*c(c1,N+3);% Feed/solution equilibrium

%Feed solution/membrane mass transfer coefficient for each component------(22)%

(Mc(c1)-J_v)*c(c1,2)+J_v*c(c1,N+3)+(z(c1)*c(c1,2)*Dinf(c1)*(f/(gc*T)))*zeta-Mc(c1)*c(c1,1);%Thermodynamic equilibrium conditions at feed solution/membrane for each component--------(24)

(gamma(c1,3)*c(c1,3))+(gamma(c1,2)*c(c1,2)*z(c1)*(f/(gc*T))*PHI(c1)*PHI_B(c1)*exp((-z(c1)*f*PSI(3))/(gc*T))*PSI(3))-(gamma(c1,2)*c(c1,2)*PHI(c1)*PHI_B(c1)*exp((-z(c1)*f*PSI(3))/(gc*T))*(1+(((z(c1)*f*PSI(3))/(gc*T)))));%Thermodynamic equilibrium conditions at membrane/permeate-solution interface for each component----(33)

(gamma(c1,N+2)*c(c1,N+2))-(PHI(c1)*exp(z(c1)*(f/(gc*T))*(PSI(N+3)-PSI(N+2))))*c(c1,N+3)+(gamma(c1,N+3)*c(c1,N+3)*z(c1)*(f/(gc*T))*PHI(c1)*exp(z(c1)*(f/(gc*T))*(PSI(N+3)-PSI(N+2))))*PSI(N+2)-(gamma(c1,N+3)*c(c1,N+3)*z(c1)*(f/(gc*T))*PHI(c1)*exp(z(c1)*(f/(gc*T))*(PSI(N+3)-PSI(N+2))))*PSI(N+3)-(gamma(c1,N+3)*c(c1,N+3)*PHI(c1)*exp((z(c1)*(f/(gc*T)))*(PSI(N+3)-PSI(N+2))))*(z(c1)*(f/(gc*T))*(PSI(N+3)-PSI(N+2)));%Electroneutrality at each node--------(32)"

z(1)*c(1,c2)+z(2)*c(2,c2)+c_x;%Electroneutrality in boundary layer on feed solution/membrane------(23)

z(1)*c(1,2)+z(2)*c(2,2)+z(3)*c(3,2);%Electroneutrality in boundary layer on permeate solution/membrane------(34)

z(1)*c(1,N+3)+z(2)*c(2,N+3)+z(3)*c(3,N+3);];endendend

## Best Answer