% Find the solutions to x^2 +y^2 +z^2 =4\;\;\;\;x+\mathrm{x6}+z=1\;\;\;\;x-y-z=1
%

syms x y z
eq1 = x^2 + y^2 + z^2 == 4;
eq2 = x + x*y + z == 1; 
eq3 = x - y - z == 1;
[xsol, ysol, zsol] = solve([eq1,eq2,eq3],[x,y,z])

%
% Plot P\left(x,\sigma \right)=\frac{1}{\pi x}\left\lbrack \frac{\sigma }{{\mathrm{log}\left(x\right)}^2 +\sigma^2 }\right\rbrack for \sigma =0\ldotp 1,\;\;\;\sigma =0\ldotp 5,\;\;\sigma =1  and 0\ldotp 01<x<3
%

syms x sigma
syms P(x,sigma)
P(x,sigma) = sigma/(log(x)^2+sigma^2)/(pi*x)
figure
fplot(P(x,0.1),[0.01,3])
hold on
fplot(P(x,0.5),[0.01,3])
fplot(P(x,1),[0.01,3])
axes1 = gca;
% Create xlabel
xlabel('x');
box(axes1,'on');
axis(axes1,'tight');
% Set the remaining axes properties
set(axes1,'FontSize',14);

%
% Find \int_0^x P\left(\xi ,\sigma \right)d\xi        (assume \sigma is real)
%

assume(sigma,'Real')
syms xi
int(P(xi,sigma),xi,[0,x])

%
% Find the value of x that maximizes P\left(x,\sigma =0\ldotp 5\right)
%

eq4 = diff(P(x,0.5),x) ==0;
xatstationarypoints = solve(eq4,x)
xatmax = xatstationarypoints(2)

%
% Find the maximum value of P\left(x,\sigma =0\ldotp 5\right)
%

P(xatstationarypoints,0.5)
simplify(P(xatmax,0.5))
double(P(xatmax,0.5))

%
% Plot P\left(x,\sigma \right)=\frac{1}{\pi x}\left\lbrack \frac{\sigma }{{\mathrm{log}\left(x\right)}^2 +\sigma^2 }\right\rbrack for \sigma =0\ldotp 1,\;\;\;\sigma =0\ldotp 5,\;\;\sigma =1  and
%

clear all
syms x sigma
P = sigma/(log(x)^2 + sigma^2)/(pi*x)
Pfor01 = subs(P,sigma,0.1)
figure
fplot(Pfor01,[0.01,3])
hold on
Pfor05 = subs(P,sigma,0.5)
fplot(Pfor05,[0.01,3])