I thing there are two problems with your approach. One is that you could simplify your code like this:
WL= [ 6.79 6.89 6.38 5.67 6.91 7.66 6.41 6.04 6.41 5.55 6.93 5.67 6.69 6.51 6.32 7.33 6.43 6.52 6.13 6.72 6.7 7.78 6.18 5.41 6.94 6.51 6.02 5.94 6.08 ];
% observed water level data
pd = fitdist(WL','ev');
% Creates a ProbDistUnivParam object by fitting the data to a extreme value distribution.
ci = paramci(pd);
% This function calculates the values of the parameters based on a certain confidence interval. Here the by default the confidence interval is 95 percent
probplot('extreme value' , WL);
h=probplot(gca,@evcdf,[ci(1,1),ci(1,2)]);
set(h,'color','r','linestyle','-')
t= probplot(gca,@evcdf,[ci(2,1),ci(2,2)]);
set(t,'color','g','linestyle','-.')
Note that probplot is by default a normal distribution. You should specify whenever you use a different distribution. Secondly, what makes you think that there is an error with what you did? This is not curve fitting, you are comparing cdfs. There is no reason why the parameters you obtained from ci will give you two cdf's "above" and "below" the best fit cdf. You probably expect that from vanilla curve fitting. But these are not the "curves" you fitted to, these are cdf's.
Please take a look at the cdf's, as they might help you understand better:
WL= [ 6.79 6.89 6.38 5.67 6.91 7.66 6.41 6.04 6.41 5.55 6.93 5.67 6.69 6.51 6.32 7.33 6.43 6.52 6.13 6.72 6.7 7.78 6.18 5.41 6.94 6.51 6.02 5.94 6.08 ];
% observed water level data
pd = fitdist(WL','ev');
% Creates a ProbDistUnivParam object by fitting the data to a extreme value distribution.
ci = paramci(pd);
% This function calculates the values of the parameters based on a certain confidence interval. Here the by default the confidence interval is 95 percent
evBest = @(x) evcdf(x,pd.Params(1),pd.Params(2));
evLow = @(x) evcdf(x,ci(1,1),ci(1,2));
evHigh = @(x) evcdf(x,ci(2,1),ci(2,2));
fplot(evBest,[0 10],'k--');
hold on
fplot(evLow,[0 10],'r-');
fplot(evHigh,[0 10],'b-');
