solve Error using mvnrnd "SIGMA must be a symmetric positive semi-definite matrix" in line :"wn = mvnrnd(mN,SN)"

%apply Bayesian Linear Regression to reach 3.8 & 3.9 pics of...
%bishop book on data which X is input and T is target...
%vector(gold standard, desire output) with number N2.
%% Initialization
clear 
close all
clc

N = 100;				% training length----100-OK
M = 10;				% model order limit-ok-0.1*100  %M is total number of parameters wj-OK
alpha 	= 2;          %p.154-OK
beta	= 25;         %p.154-OK
lambda  = alpha/beta;		%p.153-OK -regularization factor
sgm = .1;			% Gaussian kernel parameter-OK
%%  a data set, a sine function data 
load('DataSet.mat');x=X;
uj 	= linspace(0,1,M)'; 				% M in number of basis function[10*1]
y   =T;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 9 Gaussian basis functions%%%%%%%%%%%%%%%%%
X = repmat(x,1,length(uj));%[100*10]
U 		= repmat(uj',length(x),1);%[100*10]=[length(x),length(uj)]
Phi_x	= exp(-(X-U).^2/((2*sgm).^2));	% (3.4)-basis function corresponding to x
wn 	  	= normalRegEqn(Phi_x,y,lambda);	% (3.28)-p.145-basis coefficients-weghts
%%
a = (randperm(numel(x),25))';  % select random--> a = randi([1  100], 25 , 1);
x1=x(a);    % matrix randomly
y1=y(a);

X 		= repmat(x1,1,length(uj));%[25*10]
U 		= repmat(uj',25,1);%[25*10]
phi_xs	= exp((-(X-U).^2)/((2*sgm).^2));		%[25*10]% basis corresponding to x


SN 	   	= inv(alpha*eye(M)+beta*(phi_xs)'*phi_xs);%[10*10]
mN	   	= beta*SN*phi_xs'*y1;%[10*1]

W = [];
for lp = 1:100;
% 	wn		= mvnrnd(zeros(M,1),eye(M)/alpha)';
    wn		= mvnrnd(mN,SN);%%%%??????????????????????????????????????????????????????????????????????????????????????????????
	W 		= [W wn];
end

plotY(x, Phi_x, W, uj, sgm);
plot(x1,y1,'bo','LineWidth',3,'MarkerSize',16);
%%