clc;
clear;
Pr=1;
R0plus = input('Enter a value for R0+: ');
syms yplus Yplus Yp;
rplus = R0plus-yplus;
Eoverv = piecewise(yplus<=50, (0.4.*(yplus-(11*tanh(yplus./11)))), yplus>50, (0.4/3)*(yplus.*(0.5+(((R0plus-yplus)./R0plus).^2))).*(1+((R0plus-yplus)./R0plus)) );
up = (R0plus-yplus)./(Eoverv-1);
q = int(up,yplus, 0, Yplus);
uplus = (q./R0plus);
um = (uplus.*(R0plus-Yplus))/(R0plus^2);
Umplus = vpaintegral(um, Yplus, 0, R0plus)
ReD= 2*Umplus*R0plus;
syms ypprime yplusprime
third = @ (ypprime) (uplus.*(R0plus - ypprime));
three = int(third,ypprime, [yplusprime, R0plus]);
secnd = @(yplusprime) three./((R0plus-yplusprime).*((Eoverv)+(1/Pr)));
b= int(secnd,yplusprime, [0 Yplus]);
first= @(Yplus) ((R0plus-Yplus).*b);
c= int(first,Yplus, [0 R0plus])