Calculate the variance of $f_X(x)$

We have $$f\left(x\right)=0.8\left(\frac{1}{5\sqrt{2\pi}}e^{-\frac{1}{50}\left(x-50\right)^2}\right)+0.2\left(\frac{1}{8\sqrt{2\pi}}e^{-\frac{1}{128}\left(x-60\right)^2}\right)$$

This can be written in a more nicer way:

$$g\left(x\right)=0.8\left(\frac{1}{5\sqrt{2\pi}}e^{-\frac{1}{2}\left(\frac{x-50}{5}\right)^2}\right)+0.2\left(\frac{1}{8\sqrt{2\pi}}e^{-\frac{1}{2}\left(\frac{x-60}{8}\right)^2}\right)$$

To calculate the expectation, it is equal to $$0.8(50) + 0.2(60)=52$$, since we see that this is essentially $$0.8N(\mu=50, \sigma^2=25)+0.2N(\mu=60, \sigma^2=64)$$

We get confirmation from desmos:

$$var(0.8N(\mu=50, \sigma^2=25)+0.2N(\mu=60, \sigma^2=64)=0.8^2(25)+0.2^2(64)=18.56$$

But also, $$v(x)=E(x^2)-(E(x))^2$$

Verifying the above on desmos, I get a different answer:

About $$48.8$$, so I am not sure where I went wrong. I know the integral is entered correctly on desmos, but where have I gone wrong with the variance?

$$f_X$$ is the density of the random variable $$UZ_1+(1-U)Z_2$$ where $$U\sim Bernoulli(0.8), Z_1\sim N(50,25)$$ and $$Z_2\sim N(60,64)$$ and they are mutually independent.

Then $$Var(UZ_1+(1-U)Z_2)=E[Var(UZ_1+(1-U)Z_2|U)]+Var[E(UZ_1+(1-U)Z_2|U)].$$

$$E[Var(UZ_1+(1-U)Z_2|U)]=E[U^2Var(Z_1|U)+(1-U)^2Var(Z_2|U)]=E[U^2]Var(Z_1)+E[(1-U)^2]Var(Z_2)=E[U]\times 25+E[1-U]\times 64=0.8\times 25+0.2\times 64=32.8.$$

$$Var[E(UZ_1+(1-U)Z_2|U)]=Var[UE(Z_1)+(1-U)E(Z_2)]=Var[50U+60(1-U)]=Var[-10U]=100\times 0.8\times 0.2=16.$$

for a function of the form

$$f(x) = a \frac{1}{\sqrt{2\pi \sigma_1^2}}e^{(x - \mu_1)^2/2\sigma_1^2} + (1- a)\frac{1}{\sqrt{2\pi \sigma_2^2}}e^{(x - \mu_2)^2/2\sigma_2^2}$$

you have

$$\int_{-\infty}^{+\infty} {\rm d}x~ f(x) = a + (1 - a) = 1$$

$$\int_{-\infty}^{+\infty} {\rm d}x~ x f(x) = a\mu_1 + (1 - a)\mu_2 \tag{1}$$

and

$$\int_{-\infty}^{+\infty} {\rm d}x~ x^2 f(x) = a(\mu_1^2 + \sigma_1^2) + (1 - a)(\mu_2^2 + \sigma_2^2) \tag{2}$$

So that

$$\mathbb{V}{\rm ar}[X] = \int_{-\infty}^{+\infty} {\rm d}x~ x^2 f(x) - \left( \int_{-\infty}^{+\infty} {\rm d}x~ x f(x) \right)^2 = a(\sigma_1 - \sigma_2)(\sigma_1 + \sigma_2) + \sigma_2^2 + a(1 - a)(\mu_1 - \mu_2)^2 \tag{3}$$