How to take the integral of a symbolic function (with/after substitutions) and plot the result?

So what I'm trying to do is solve a symbolic ODE, then use the results (after making some substitutions) to compute various things and plot them. So far everything works fine, but what I'm trying to do now is take the integral of one of my symbolic functions (power) to get energy, and plot the result vs. time, which is not going so well.

What I've tried so far is this:

1. Taking the integral of my power function using int(), then plotting the result with fplot like this:

%taking the integral:
energyinmJ = (int(powerinW,t))*1000;
%and the line for the plot:
fplot(energyinmJ,'--o','LineWidth',2)

This results in the computation running a very long time (it normally takes about 30 seconds, this never finished and I let it run for about 30 mins). I also tried taking the integral over just [0,20], which (this is perhaps a clue) results in many lines instead of just the one line I am expecting.

2. I also tried converting the power function to a vector of values over the interval from 0-20, then using cumtrapz to approximate the integral, which now almost works, but not quite (see below).

I am obviously doing something wrong, but I am at the point where I could use some help. What I am expecting for the integral is a single stair-step-like function that starts at zero and increases with time (and proportionally to "subsdiffsolnthetadot"). Any help would be appreciated. If there's just a better method to use and I need to start over, I'm open to that, too. Full code posted below. Thanks!

-J

Edit: I am closer to a result using cumtrapz, by doing this:

FUN = matlabFunction(powerinW);  
y = feval(FUN, t);
energyinJ = cumtrapz(t,y);

But despite t being defined by t = 0:0.1:20; energyinJ ends up being an array 201 units long, with only 20 actual values and the rest NaN. Continuing to look for the problem, but still, any insight would be appreciated. Thanks again! Edited code posted below.

-J

     %Set up all variables used in symbolic ODE:
  syms theta(t) m r EGR Q u k taustall A B C t; 
  thetadot(t) = diff(theta(t),t);

%Input symbolic ODE, initial conditions, solve with dsolve:
ode = 0.5*m*(r^2)*diff(theta(t),t,2)+(((-Q*EGR)+u)*diff(theta(t),t))+(k*(r^2)*theta(t))+(taustall*EGR) == (((A/2)*sin((B*t)+C))+(A/2))*r
pretty(ode)
cond1 = theta(0) == 0;
cond2 = thetadot(0) == 0;
conds = [cond1 cond2];
solntheta(t) = dsolve(ode,conds);

%Simplify solution for theta and take first derivative to obtain angular velocity:
solntheta = simplify(solntheta);
diffsolntheta = diff(solntheta,t);

%Set up Forcing Function as standalone for substitution & plotting:
forcingfunc = (((A/2)*sin((B*t)+C))+(A/2));

%Apply values to all variables in ODE:
m = 0.1; %about 100g for the mass of the pulley
r = 0.025; %meters (about 5cm in diameter)
EGR = 1; %because no external gearing yet, only whatever's in the gearmotor (and therefore accounted for by the torque-speed curve)
Q = (0.003531/513.3); %stall torque in N-m over rated rpm in rad/sec 
u = 0.01; %assume very low friction
k = 30.7; %30.7 N/m from that paper
taustall = 0.003531; %N-m from motor specs sheet, converted from 0.5 oz-in
A = 0.45; %assume minimal force of chest expansion from paper (Max amplitude is actually A/2 for this func)
B = 6/(2*pi); % Because T=2pi/B, therefore B=T/2*pi, Bslow=6/2*pi = 0.1Hz which is 10 breaths/min (T=6sec), or use Bfast = 2/2*pi = 0.5Hz which is 30 breaths/min (T=2sec)
C = (3*pi)/2; %For a phase shift of 0, (3*pi)/2 when B=T/2*pi.

%Set desired interval to calculate over here:
%t = linspace(0,20,300); %plot over 20 seconds (about 3 breath cycles)
t = 0:0.1:20;

%Homogenous side terms (for reference only)
%a = 0.5*m*r^2;
%b = ((-Q*EGR)+u);
%c = k*r^2;
%d = taustall*EGR;

%Make substitutions and store results:
subsofsolntheta = subs(solntheta);
subsofdiffsolntheta = subs(diffsolntheta);
subsofforcingfunc = subs(forcingfunc);

%POWER CALCULATIONS:
%Power = torque*angularvelocity = I^2*R
powerinmW = abs(1000*(((-0.003531/513.3).*subsofdiffsolntheta)+0.003531).*subsofdiffsolntheta);
powerinW = abs((((-0.003531/513.3).*subsofdiffsolntheta)+0.003531).*subsofdiffsolntheta);
%P=IV=I^2*R b/c V=IR so...
%I=sqrt(P/R), set R so I<=1600mA and V high enough to use, in this case 3000 Ohm works for now.
currentinmilliAmps = 1000*(sqrt(powerinW./3000)); 
%And V=IR, with R same as above.
Voltage = (currentinmilliAmps/1000)*3000;

%Create function handle for power then use cumtrapz():
FUN = matlabFunction(powerinW);  % This creates a function handle
y = feval(FUN, t);
energyinJ = cumtrapz(t,y);

%PLOT EVERYTHING:
figure
hold on
plot(t,subsofforcingfunc,':k','LineWidth',3)
plot(t,subsofsolntheta,'-b','LineWidth',2)
plot(t,subsofdiffsolntheta,'-m','LineWidth',2)
plot(t,powerinmW,'--r','LineWidth',2,'MarkerIndices',1:15:length(subsofsolntheta))
plot(t,currentinmilliAmps,'--*c','LineWidth',3,'MarkerIndices',1:15:length(subsofsolntheta))
plot(t,Voltage,'--gv','LineWidth',3,'MarkerIndices',1:15:length(subsofsolntheta))
title('Forcing function, theta(t), thetadot(t), Power, Current, & Voltage (rectified) vs. t');
xlim([0 20]);
xlabel('t'); 
ylabel({'Forcing function, theta(t) (in radians), thetadot(t) (in radians/sec)';'Power (mW), Current (mA), Voltage (V)'});
legend('Forcingfunc','theta(t)','thetadot(t)','Power in mW','Current in mA','Voltage in V','Location','northeast');
hold off

figure
hold on
plot(t,powerinmW,'--r','LineWidth',2,'MarkerIndices',1:15:length(subsofsolntheta))
plot(t,energyinJ)
% this almost worked: syms t
% this almost worked: fplot(energyinmJ,[0,20])
title('Power in mW and Energy in J vs. t');
xlim([0 20]);
xlabel('t'); 
ylabel({'Power in mW, Energy in mJ'});
legend('Power in mW','Energy in mJ','Location','northeast');
hold off

Answer 1

[JustAnotherHumanoid]

Alright, I got it to work, so I'm answering my own question here in case someone else has a similar problem. I started over, and defined my power equation like this, after substitutions for everything except "t":

alternatepowerinW = abs((((-0.003531/513.3).*subs(diffsolntheta))+0.003531).*subs(diffsolntheta));

Then I set up the following variables, which I named to deconflict with other things I was trying out at the time:

powerfunc = alternatepowerinW;
pointstoeval = 0:0.1:20;
powerarray = powerfunc(pointstoeval);
energyinmJ = cumtrapz(pointstoeval,powerarray)*1000;

This took a long time (~5mins) to compute, but it produced what I expected to see, and I get actual values instead of many NaN elements (still need to verify the result, which I attached as an image). If you looked at my post and tried to help out, thanks, I'm sure it won't be long before I'm back with another question.

-J