MATLAB: I am having problem with the code written below. It is a code for finding stored strain in a polymer. Some of the errors that are coming are index exceeds matrix dimension, matrix dimensions must agree, etc. Kindly help me out.

polymer

syms epre et es(T); %stored strain
Tl = 273;
Th = 358;
Tg = 343;
Ei = 813E6; %modulus corresponding to the internal energetic deformation
N = 5.938E26; 
K = 1.38E-23; %[J/K] Boltzmann's constant
a = 2.76E-5;
n = 4;
es = zeros(1,100);
epre = .091; %pre-strain under tension
Phif = 1-1./(1+a*(Th-T).^n);
y = diff(Phif,T);
fun = @(T)1.42E-6.*T-3.16E-4;
for T = Tl:Th 
   
    Ee = 3*N*K*T;
    E = 1./((Phif./Ei)+((1-Phif)./Ee));
    et = integral(fun,Th,T); 
    
    Des = diff(es);    
    ode = Des == (epre-es-et)*y.*E./Ee;%writing stored strain as O.D.E
    cond = es(Th) == 0; %initial condition for O.D.E
    es = dsolve(ode,cond); %solving O.D.E
    
    disp(es);
    
    plot(T/343,es,'.'); hold on
end

Best Answer

Try this: 
syms epre et es(T); %stored strain
Tl = 273;
Th = 358;
Tg = 343;
Ei = 813E6; %modulus corresponding to the internal energetic deformation
N = 5.938E26; 
K = 1.38E-23; %[J/K] Boltzmann's constant
a = 2.76E-5;
n = 4;
epre = .091; %pre-strain under tension
    Phif = 1-1./(1+a*(Th-T).^n);
    
    y = diff(Phif,T);
    
    fun(T) = 1.42E-6.*T-3.16E-4;
   
    Ee = 3*N*K*T;
    E = 1./((Phif./Ei)+((1-Phif)./Ee));
    et = int(fun,Th,T); 
    
    Des = diff(es);    
    ode = Des == (epre-es-et)*y.*E./Ee;%writing stored strain as O.D.E
    cond = es(Th) == 0; %initial condition for O.D.E
    [VF,Sbs] = odeToVectorField(ode)
    esfcn = matlabFunction(VF)
    
    Tv = linspace(Tl, Th);
    [T,Es] = ode45(esfcn, Tv, 0);
    
    
    figure
    plot(Tv, Es)
    grid
    xlabel('T (°K)')
    ylabel('Es')
This does not give you an analytical solution (that may not actually exist), however it does give you a plot of your result.
Related Question