Gaussian Process Binary Classification: Calculating Predictive Output Variances on Probability Scale

I'm attempting to use MATLAB's gpml toolbox ( http://www.gaussianprocess.org/gpml/code/matlab/doc/ ) to do classification and I'd like to have confidence bands around the final predicted probabilities. I can get the confidence bands around the latent function and/or the predictive output means. With binary classification, however, the predictive output mean first has to be transformed to probabilities. Using MATLAB 2013a, I have the following:

%-------------------------------
% Create data from test cases
n = 30;
x = 10 * lhsdesign(n, 1);
prob_fun = @(x) 0.75 * normcdf(-x,-1.75,0.4) + 0.5 * normpdf(x,4.5,1) + 1.75*normpdf(x,7.5,0.75);
prob = prob_fun(x); 
y = binornd(1, prob, n, 1);
test_cases = linspace(min(x), max(x), 500)';

% Convert to -1/1 for gp code. Also see what true function looks like. 
y(y < 1) = -1;
true_probability = prob_fun(test_cases); 
% plot(xs, truth, 'k-')

%-------------------------------
% Set mean to be constant. Put in terms of logit
meanF = {@meanConst};
meanY = mean(0.5 * (y + 1));
meanY = log(meanY / (1 - meanY));
hyp0.mean = meanY;

% Gaussian correlation (covariance) function. Just manually setting the
% length and scale parameters for now
covfunc = @covSEiso;
hyp0.cov  = [0.75; 2.5];
likfunc = @likLogistic;

% Run GP model and make predictions on test cases
[ymus, ys2s, fmus, fs2s, ~, post] = gp(hyp0, @infEP, meanF, covfunc,...
    likfunc, x, y, test_cases);

%-------------------------------
% Turn the probability values into valid probabilities:
ymus_prob = (ymus + 1) * 0.5;    

% THIS IS WHERE I'M STUCK...IS THIS CORRECT?
ys2s_lower_prob =  normcdf(ymus + 1.96* sqrt(ys2s));
ys2s_upper_prob =  normcdf(ymus - 1.96* sqrt(ys2s));  

% Realizations converted
y_01 = y;
y_01(y_01 < 0) = 0;

%-------------------------------
% Plotting
figure()
subplot(1, 2, 1);
plot(x, y, 'ko'); hold on;                                      % realizations
f = [ymus + 1.96*sqrt(ys2s);...
flipdim(ymus - 1.96*sqrt(ys2s), 1)];
fill([test_cases; flipdim(test_cases,1)], f, [7 7 7]/8);        % confidence region
plot(x, y, 'ko'); hold on;                                      % realizations
plot(test_cases, true_probability, 'k--'); hold on;             % true function
plot(test_cases, ymus, 'r-'); hold on;                          % predicted function 
title('Predicted Values: Not Transformed')

subplot(1, 2, 2);
plot(x, y_01, 'ko'); hold on;                                   % realizations
f = [ys2s_lower_prob;...
flipdim(ys2s_upper_prob, 1)];
fill([test_cases; flipdim(test_cases,1)], f, [7 7 7]/8);        % confidence region
plot(x, y_01, 'ko'); hold on;                                   % realizations
plot(test_cases, true_probability, 'k--'); hold on;             % true function
plot(test_cases, ymus_prob, 'r-'); hold on;                     % predicted function 
title('Predicted Values: Transformed to 0-1 Scale')

You can see that I'm stuck trying to figure out what to do with the ys2s and how to get it to "probability terms", so to speak. The plots are:

Can someone offer some direction on how to generate predicted variances on the probability scale? I think I did it correctly here and while I understand that the confidence bands may not be symmetric, they're not even containing the mean in some spots.

If it makes any difference, I'm on a Windows 10 machine. I also asked this question on StackExchange last week but have yet to hear anything back ( https://stackoverflow.com/questions/51730293/gaussian-process-binary-classification-calculating-predictive-output-variances )