RMS value of non-sinosoidal current

I'm trying to get an RMS value of a non-sinusoidal current function and cannot make it work in MATLAB. Here is the formula I would like to implement via MATLAB in my code. Please help implement formula (1) for obtaining the RMS value of non-sinusoidal current. Another approximate variant is the formula (2) with the summation symbol:

U = [30,60,90,120];
S = 10*10^(-4);
Phim = U/(4.44*50*400);
Bm = Phim/S;
wt = 0:3.14/50:2*3.14;
%Givens 
R = 0.5/400;
H1 = 200;
B1 = 1.0;
H2 = 1000;
B2 = 1.43;
%H1=a*B1+b*(B1)^5
%H2=a*B2+b*(B2)^5
%Solve 
syms a b 
subs(a*B1+b*(B1)^5, B1, 1)
subs(a*B2+b*(B2)^5, B2, 1.43)
[a,b] = solve(a*1+b*(1)^5 == 200, a*1.43+b*(1.43)^5 == 1000);
A = sin(wt');
B = Bm.*sin(wt');
i = ((a*B)+(b*(B.^5))).*(R);
%Plot
subplot(121)
plot(wt, B)
xlabel('\omegat'), ylabel('B, Magnetic density')
title('Магнитная Density for [30, 60, 90, 120V]')
subplot(122)
plot(wt, i)
xlabel('\omegat'), ylabel('i, А')
title('Non-sinusoidal current for [30, 60, 90, 120V]')

%Attempt at RMS value of current 
f=inline('0.811*sin(wt)-0.367*sin(wt)+0.073*sin(5*wt)');
a=0;
b=3.14/2;
h=0.0002;
N=((b-a)/h)+1;
i=1:N; 
x=a:h:b; 
y=sqrt(1/6*(feval(f,x)).^2); 
m=2:N;
y1(m-1)=y(m);
Fr=sum(h*y1);

Answer 1

[dpb]

r=rms(x);

[Krivar Afron]

Thanks for the answer. I tried it like that, but I did not exactly get the answer in the book. The answer in the book for 30V, 60V, 90V, and 127V are the values: 0.03, 0.071, 0.17, and 0.57A, respectively. What I get in MATLAB are the values: 0.0131A, 0.0380A, 0.1382A, 0.6251A, which is different than the book answer. The code I used is below. Basically, the idea is to get an RMS value of these non-sinusoidal current waveforms, each representing a certain input voltage value, from 30V to 127V. Any idea why this is happening?

r = rms(i);
double(r);
p = rms(0.811*sin(wt)-0.367*sin(3*wt)+0.073*sin(5*wt));