Problems with pdepe and input variables

I have played around with my code and have the following code for pdepe which looks as if it should work:

function battery_model

%This is the simplistic model of a lithium ion battery model which I
%cobbled together knowing very little about batteries. The constants and so
%forth come from experiment.
%Define some global constants and parameters:
global rho; rho=1;global c_th; c_th=1;
global a_1; a_1=1;global a_2; a_2=1;global b; b=1; global c_1; c_1=1;global c_2; c_2=1;
global I_app; I_app=1.5;
global k;global D; global epsilon;
k=1; %Thermal conductivity, most likely a function.
epsilon=1; %Electrical conductivity, most likely a function
D=1; %lithium ion diffusion, most likely a function
m=0;
x=linspace(0,1,100);
t=linspace(0,1,50);
%Initial conditions;
global T_0; global c_0; global phi_0;
T_0=ones(1,100);
c_0=smooth_step(x,0,1/3);
E_0=I_app/(epsilon)-(b/epsilon)*cumtrapz(x,c_0);
phi_0=cumtrapz(x,E_0);

sol = pdepe(m,@battery_GE,@battery_ic,@battery_bc,x,t);
T = sol(:,:,1);
E = sol(:,:,2);
c = sol(:,:,3);

figure;
surf(x,t,T);
title('Temperature');
xlabel('Distance x');
ylabel('Time t');

figure;
surf(x,t,E);
title('Electric Field');
xlabel('Distance x');
ylabel('Time t');

figure;
surf(x,t,c);
title('Li ion concentration');
xlabel('Distance x');
ylabel('Time t');

function [CC, FF, SS]=battery_GE(x,t,u,DuDx)

CC=[rho*c_th; 0; 1];
FF=[a_1*k; epsilon; c_1*D].*DuDx;
SS=[a_2*sigma*DuDx(2); -b*u(3);c_2*(DuDx(3)-(b*u(3)/epsilon)-DuDx(1))];

function u0=battery_ic(x)
global T_0; global c_0; global phi_0
u0=[T_0;phi_0;c_0];

function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
global I_app; global sigma; global epsilon;
T_L=1;
T_R=2;
pl=[ul(1)-T_L; ul(2); 0];
ql=[0; 0; 1];
pr=[ur(1)-T_R; -I_app; 0];
qr=[0; sigma/epsilon; 1];

I get the following error which is baffling to me:

Error using pdepe (line 231)
Invalid output of ICFUN. ICFUN must return a column vector.

Error in battery_model (line 24)
sol = pdepe(m,@battery_GE,@battery_ic,@battery_bc,x,t);

Any idea what on earth is going wrong?

[Torsten]

T_0 is a 1x100 row vector. If you want to supply it in battery_ic as you do, it should be a scalar.

Best wishes

Torsten.

[Matthew Hunt]

In general it should be a function of position, so not a scalar.

[Torsten]

"battery_ic" should return the values of the solution at time t=0 at a fixed position x, thus a 3x1 vector in your case.

[Matthew Hunt]

Is there a standard way to do this?

[Torsten]

If

X=linspace(0,1,100)

from "function battery_model", then

function u0=battery_ic(x)
global X; global T_0; global c_0; global phi_0
T_00 = interp1(X,T_0,x)
phi_00 = interp1(X,phi_0,x)
c_00 = interp1(X,c_0,x)
u0=[T_00;phi_00;c_00];
end

[Matthew Hunt]

It now works completely.

The initial conditions seems to have a bit of a faff to do. Is this the general thing which I should do?